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Preface 


‘No modern applied mathematician, physical scientist, or engineer can be properly trained without some 
understanding of numerical methods’ 
—E. Issacson and H. B. Keller 


Computational hydraulics is one of the many fields of science for which computers opened a new way 
of working. The user of computer based models should be able to critically evaluate results given by 
different computational methods used in solving practical hydraulic problems. As such the book is 
intended to serve as a course for both undergraduate and graduate students in the fields of hydraulics and 
water resources. 

The text of this book is based on the course of “Numerical methods for free surface flow", which 
evolved over the years at UNESCO-IHE Institute for Water Education. My students have been a source of 
motivation in trying to explain and clarify complex aspects of numerical methods, as well as free surface 
flows and hydrodynamics of large water bodies such as lakes and reservoirs. 

There is a great deal of literature devoted to computational hydraulics in terms of numerical methods 
and fluid flow phenomena. However, it would be impossible to cover everything in one single book. 
Hence I have tried to cover only those methods that are most used in hydrodynamics. In Chapter 1, the 
book introduces the concept of modeling and the contribution of both numerical methods and numerical 
analysis to modeling. A description of the basic hydraulic principles, and the problems addressed by 
these principles in the aquatic environment is presented in Chapter 2, followed by the discretization 
principles of the fluid flow domain in Chapter 3. The finite difference method and finite volume method 
are presented in Chapters 4 and 5, respectively, with all the necessary steps for building and applying 
these numerical method approaches in hydraulics. Chapter 6 presents the properties of numerical 
methods that gives insight into how to avoid wrong results. The final two chapters show different 
example applications of computational hydraulics: river system modeling alongside hydrodynamic and 
water quality modeling of lakes and reservoirs. Worked-out examples that help in understanding the 
main concepts of computational hydraulics are developed for demonstration purposes at the end of 
Chapters 4, 5, 6 and 7. 
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Chapter 1 
Modelling theory 


1.1 CONTEXT AND NATURE OF MODELLING 


Models come in many different forms, some of which we do not usually refer to as models. Language, for 
example is a model, making associations between abstract concepts and labels (defined by words), using a 
set of rules for clustering words together (grammar), as such allowing to describe reality (i.e., build models 
of reality). 

A model is a simplified, schematic representation of the real world. Models are meant to help engineers, 
scientists and decision-makers to determine what is happening in reality and to predict what may happen 
in the future. In particular, they are useful for the assessment of the impact of human activities on the 
environment or on artificial systems. Models are covering a large range of problems, hence it is difficult to 
find one that fits all problems. 

A classical definition of the model is ‘a simplification of reality over some time period or spatial extent, 
intended to promote understanding of the real system’ (Bellinger, 2004) or ‘A model is a simplification 
of reality that retains enough aspects of the original system to make if useful to the modeller' (Chapra & 
Canale, 2006; Eykhoff, 1974). In this context the system is defined to be a part of reality (isolated from the 
rest), which consists of entities that are in mutual relationships (processes) and have limited interactions 
with the reality outside of the system. Examples of systems are ecosystem or a hydrological system. A 
model is a physical or mathematical description of a physical system including the interaction with the 
outside surrounding environment, which can be used to simulate the effect of changes in the system itself 
or the effect of changes due to conditions imposed on the system. 

Modelling implies the construction of a model, or working with an existing model, while simulation 
means to mimic the system on a computer. Building a small prototype of a system, in order to test 
theories and design techniques was one of the first uses of a model in civil engineering and it is referred 
to as physical modelling. As science and technology progressed, so did the knowledge about systems and 
properties of a system could be quantified by mathematical equations that were either empirical or derived 
from basic principles. The exact solution of such equations is referred to as analytical model. Mathematical 
representations in form of differential equations do not always have explicit analytical solutions and the use 
of approximate numerical approaches, to solve complex equations, is referred to as numerical modelling. 
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Modelling is as much an art as it is a science. The purpose to construct a model is to gain insight into 
the underlying causes or to provide accurate prediction of the system behaviour. Reliability on the model 
results depends on the expertise of the modeller about the system to be modelled. In general models are, 
as mentioned, mathematical representations, in forms of equations, of physical phenomena, embedding 
parameters and quantitative relationships between different variables and parameters. They should be seen 
as condensed versions of the knowledge of the modeller about the modelled system. 

Approaches to construct a model are different from one problem to another as well as from one 
modeller to another. An engineer approach to problem solution is different from that of a mathematician. 
Mathematicians are interested in finding whether a solution to a differential equation exists, while an 
engineer assumes that the existence of a physical system is proof enough for the existence of a solution and 
focuses in finding the solution itself. The understanding of an engineering system is essentially gained by 
Observation and experiment. After years of observation and measurements, engineers notice that certain 
aspects of their studies occur repeatedly. Such general behaviour can then be expressed as fundamental 
laws or models. These kinds of models embody the cumulative wisdom of past experience. As the purpose 
of modelling is to increase the understanding of a system behaviour, the validity of it does not reduce to 
its fit to observations and measurements, but also to model’s ability to simulate and replicate situations or 
data beyond those originally described in the model (James & Burges 1982; Dooge, 1986). For example 
a descriptive model of the system is used to understand how the system could work; or to predict how an 
unexpected event could affect the system; or in some cases to try different control approaches in case a 
system is controlled by structures. 

The developer of a model considers that advances in modelling can only be attained by forming close 
links between measurement, theory and modelling (Woods eft al., 1995). Models are constantly being 
improved and advanced, which implies that users of the models should always obtain the up-to-date 
version of a model to be applied for finding solutions to an engineering problem (Schulze, 1995). 

Models in civil engineering refer to the modelling of the continuum, however the present book is 
concerned with the modelling and numerical methods applied to fluid flow problems. Computational 
fluid dynamics (CFD) is defined as the field of science that uses computers and numerical techniques 
to solve problems involving fluid flow. The term Computational Hydraulics was defined by Abbot 
and Minns (1998) to be 'the reformulation of traditional hydraulics to suit the possibilities and 
requirements of discrete, sequential and recursive processes of digital computation’, and it is a sub- 
field of CFD. 

In order to build a model apart from the mathematical formulation of the phenomena a lot of 
Observed data is required to perform simulations. For the majority of hydraulic problems there is an 
imbalance between the lack of available data (little knowledge about the problem) and the availability of 
sophisticated methods that can be applied to simulate the system under study. Numerical methods and 
simulations fill the ‘gap’ between theory and experiment, providing qualitative and quantitative insight 
into many phenomena that are too complex to be dealt with by analytical methods and too expensive, 
time consuming and/or impossible to be performed experimentally. Hence it can be stated that simple 
models provide insight, complicated models provide many results for which insight need to be carefully 
analysed. 


1.11 Classification of models 


Classification of models is done in order to understand the type and level of mathematics involved in 
developing a model and it helps understanding how the simulation is done (Lane & Nichols, 1993). 
Moreover knowing the mathematics employed the user of a model can determine the nature of the output, 
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as well as the amount of information required to build the model. Model classification helps users to select 
the appropriate model for a particular problem solving need. 

Before any classification of models is done several formal mathematical notions and definitions are 
introduced: 


* independent variables used in modelling are space and time. Time is usually defined over an interval, 

te [t,, T], and space x refers to the volume V that contains the system under study, x € V; 

* dependent variable is the state variable, which takes values depending on parameters and independent 
variables. The state variable is a finite dimensional vector u = u(x, f), which is n dimensional u = (u,, 
U>,...U,,), and describes sufficiently the evolution of the phenomena (real system); 

* amathematical model is the set of equations that defines the evolution of the state variable in space 


and time. 


In order to properly represent the phenomena in nature, the real systems and the mathematical model 
need to be consistent. This means that the number of unknown dependent variables must be equal to the 
number of independent equations. 

Purely from mathematical point of view models are dynamic or static, or they are finite and continuous. 
A mathematical model is dynamic if the state variable u is time dependent. If u does not depend on time 
then the mathematical model is static. A mathematical model is finite if the state variable does not depend 
on the space variables. Otherwise the mathematical model is continuous. Finite dynamic models are 
represented through ordinary differential equations; while continuous dynamic models are represented by 
partial differential equations. Static models, finite or continuous, are particular cases of dynamic models 
(the case when time derivative is zero). 

From engineering point of view mathematical models for water related problems can be described as: 


* Lumped conceptual models that are based on the concept of exchanges between global storage 
entities, which compose the system under study (i.e., one compartment for surface water, another one 
for aquifer storage, etc). These models satisfy principles such as the continuity principle, but do not 
embed a complete description of the driving forces, which govern the system to be modelled (e.g., 
rainfall-runoff model, reservoir model, etc). 

e Physically-based, distributed models, that use a description of the physical phenomenon which 
govern the behaviour of water in the system under study. The principles that are applied are mass 
conservation and additional laws describing the driving forces, such as momentum equation. In case 
that these models refer to flow in saturated or unsaturated porous media, along with the continuity 
principle, the hypothesis of laminar flow is applied (i.e., shear stress proportional to the velocity, 
which determines the Darcy or Richards' equations for flow). 

* Data driven models that seek for a correlation between input and output data, without trying to 
detail/analyse phenomena (e.g., linear regression, unit hydrograph). These models do need a lot of 
prior knowledge and data observations of the system under study. 


Lumped conceptual and physically-based models use differential relationships, describing the changes 
in time of a set of variables. Most of the time, the formulations used are either Ordinary Differential 
Equations (ODEs) or Partial Differential Equations (PDEs). As these ODEs and PDEs are very often 
complex, it is not always possible to find analytical solutions to them. 

Researchers (Singh, 1995; Tim, 1995), further classified models depending on the description of: 


— Problem type and solution (Figure 1.1). 
— Processes addressed (water quality, water allocation, reservoir operation, flood management, etc). 
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— Time and space problem dimension (lumped conceptual, one dimensional, two dimensional and 
three dimensional) (Figure 1.2); 
— Method of solution (finite differences, finite elements, etc. See Figure 1.3.) 


Physical system 
Problem type 
(water allocation, Water quality, river 
basin managment, flood control, etc.) 


EE DAMM M FADA MAMMA LM MM | MID NCMO: ! 


Mathematical representation of the 
physical system (Models) 


MEME 


Theoretical representation Empirical representation 


Deterministic Stochastic Hybrid 


White box | Black box Conceptual 


Grey box 


Figure 1.1 Classification of models based on model type. 


Scale 
Space Time 
| Large | 
Lumped Distributed — Medium iil term Long term 
T p mj 
ID 2D 3D Small Daily Seasonal Trends 


over years 
Figure 1.2 Classification of models based on scale. 


Mathematical modelling problems are structurally classified into black box and white box models, 
depending on how much information is available about the system under study before the formulation 
of the mathematical model to be used for the problem solution. A black-box model is a system for which 
no a priori information is available and no understanding of the processes involved in the transformation 
is required. Only the input and the output have physical meaning. Stochastic models are black box 
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models. A white-box model is a system where all necessary information is available and all processes are 
understandable and accounted for. Deterministic and physically based models fall into this category. The 
term grey box can also be used if not enough data is available and just partial understanding of processes 
is available. Lumped conceptual models are grey box models. In practice it is good to use as much a 
priori information as possible because it makes models behave correctly (accurately) (see Figure 1.4). 


Phenomena in nature 
(described in terms of properties that 
prevail at each point in time and 
space separately) 


Partial differential equations (PDEs) 
governing the phenomena 


Solution of PDEs 


—— — 


Analytical solutions Numerical techniques 
Separation of Integral Integral FDM FVM FEM 
variables solutions transforms (Finite (Finite (Finite 

(Green) (Fourier & Difference Volume Element 
Laplace) Method) Method) Method) 


Figure 1.3 Classification of models based on method of solution. 


OO LC n Knowledge on physical processes 


Mm. s "M Numerical Deterministic 

Serm TM I" models numerical 
Fullydata © 0 m ~and | models Fully process 
oriented models Hybrid Data 777 E. oriented models 


models assimilation — "7 


Measured Data 


Figure 1.4 Models coverage based on data availability. 


1.1.2 Computational Hydraulics 


The field of Hydraulics can be divided into four major areas: Theoretical, Applied, Computational 
and Experimental. Theoretical area deals with fundamental laws and principles of mechanics studied 
for their intrinsic value. The Applied area transfers the theoretical knowledge to scientific and 
engineering applications, especially as regards to the construction of mathematical models of physical 
phenomena. Computational Hydraulics solves specific problems by combining mathematical models 
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with numerical methods implemented on digital computers, a process called simulation. Experimental 
fluid mechanics puts physical laws, mathematical models and numerical simulations under the test of 
observation. Computational Hydraulics is interdisciplinary and has four major contributing disciplines 
(see Figure 1.5). 


Computational Hydraulics 


Numerical 


Hydraulics, 
methods 


Hydrodynamics, 
Hydrology 


Applied mathematics 
& analysis 


Computer & 
Information sciences 


Figure 1.5 Computational Hydraulics field of study. 


Engineering problems related to Hydraulics are of three kind; steady state, propagation problems and 
eigen-value problems. Steady state problems have as main characteristic the fact that the system does not 
change in time. Thus the equations representing the problem do not involve time as variable. Propagation 
problems are also called dynamic problems and the main characteristics of it is that the system changes 
with time. The variables and equilibrium relations depend in this case on time. The objective of the analysis 
is to calculate the state variables in time. In the case of eigen-values problems there is no unique solution to 
the system equilibrium equations. The analysis of the system has as the objective to determine the various 
possible solutions. Eigen-value problems arise in both steady-state and dynamic analysis. 


1.2 CONCEPTUALIATION: BUILDING A MODEL 


In general, mathematical modelling consists of several steps, as shown in Figure 1.6. 
The starting point for any mathematical model is the Problem Description when the definition of what 
are the physical elements that are relevant in describing the phenomenon and the processes involved are 


defined. This stage of modelling requires input from experts in computational hydraulics and sometimes 
from other disciplines. 
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Possible errors. due to modelling solution 


Physical system Mathematical formulation 
(Problem description) (Mathematical model) 
A 


Discretization via 
SF(FDM) 
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WE(FEM) 


discretization 


Discrete model 
Visualization< & p 
Numerical solution 


Possible errors due to Possible errors due to 


solution implementation 


Figure 1.6 The cycle of mathematical modelling. 


The next step is the Mathematical Formulation of the problem. After determining all relevant aspects, 
these are translated into mathematical relations that are in the form of a system of partial differential 
equations for various field quantities, which have to be solved over a given suitable domain (e.g., the whole 
volume of a lake, the whole river length). Boundary conditions are required over the domain of computation. 

The third step is the one of translating the problem from the language of mathematics into a computable 
form. This is done by generating a grid of the computational domain on which the mathematical equations 
are formulated in a discrete form using different numerical approaches such as finite differences (FDM), 
finite volume (FVM), finite elements (FEM), boundary elements (BEM), and so on. The grid generation 
itself poses challenges for complex geometries applications. The discretization method is a matter of choice 
of the modeler and in majority of the cases depends on experience, but it should correspond to the grid used. 

At the end of the discretization step a system of nonlinear or linear algebraic equations is obtained and 
appropriate solvers have to be used in order to find the solution. In case that nonlinear system of equations are 
obtained iterative methods based on successive linearization are used. This approach may lead to large systems 
of linear equations that are again solved using iterative methods. This step is known as Numerical Solution. 

After obtaining a solution, visualization of results is needed. Post-processing is used to visualize the 
simulation results. This step enables interpretation of results, which is very important not only for the 
practitioner, but also for the developer of a model, because it helps to determine if there are errors in 
the computations and to see what is the quality of the model. 

If results of simulation differ from the real situation, the model and/ or the numerical methods have to 
be refined and the cycle of mathematical modelling starts again. 


1.3 MATHEMATICAL MODELLING IN PRACTICE 
1.3.1 Selecting a proper model 


Because of the development of numerical approaches many simulation techniques have been developped 
which lead to the proliferation of different models. As already mentioned models are intended for specific 
applications, however there are many models that are very general and have a wide area of applicability. 
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Naturally the question arise which model to choose from the wide variety of available models. The choice 
is difficult because there are models that are licence free as well as commercial ones. 

Several authors (Woolhiser & Brakensiek, 1982; Sorooshian, 1991; Burnett & Watson, 1995) are defining 
practical criteria to be considered in selecting a model; the ability of the model to address the nature of the 
problem to be solved; availability of data for a particular problem; cost of the selected tool; possibility of 
applying the model for other similar problems. 

The first criteria to consider is the ability of the model to address the problem that needs solution, by 
looking at the processes that are represented in the model, what are the equations representing these and 
what are the assumptions and limitations of the mathematical model. The availability of data is a very 
important issue, because even if the model is capable to describe processes that would be interesting 
for the modeller, if there is no data available to instantiate the model then a simpler model has to be 
considered. 

Models range from very simple formula to complex ones that require implementation on a computer. 
there are advantages and disadvantages in using a complex or a simple model, but the most important is the 
objective for which the model is build. The choice to choose a complex model is an important decision in 
the modelling process, in the light data requirements, model applicability and computer power (Bergstróm, 
1991). Complex models should be used for complex simulations, such as water quantity and quality for an 
entire catchment, however they should be used with care because they require a lot of input and there is a 
danger to introduce errors (Hughes & Beater, 1987). Complex models need long computational times and 
a thorough understanding of the model. 


1.3.2 Testing a model 


Any mathematical model is applied to solve engineering problems needs to be tested through validation, 
calibration and verification. 

Schultze (1995) mentions that verification and validation notions are two terms used interchangeably by 
modellers, ‘what is one modeller's verification may be another modeller's validation’. The use of the two 
terms might be misleading (Bredenhoeft & Koniko, 1993), hence the definition of the notions, according 
to the Oxford English Dictionary is given below in order to clarify terminology: 


— Validate is defined as ‘well founded and applicable, sound and to the point, against which no 
objection can be fairly brought’ (i.e., authentic, true); 

— Verify is defined as ‘to test the accuracy, or establish the truth or correctness, of something by 
examination or by comparison with known data or some standard' (i.e., evidence, approval). 


Any mathematical model selected to seek the solution for a given problem is first validated by using it 
to simulate small scale simple cases for which the results are known or can be easily obtained, analytically 
or by measurements. The validation ensures the applicability of the mathematical model. After selection 
of the model, the particular problem case is modeled by instantiating the model with data. The evaluation 
of the performance of it will take place in two steps; calibration and verification. During the calibration 
step all model parameters are estimated based on comparison of the simulation results with observed 
data. In the calibration phase the ability of the model to reproduce the response of the system is tested, 
as it is for example the simulation of a downstream hydrograph based on an inflow hydrograph upstream. 
The next step is verification when a simulation is performed, using parameters values obtained during 
calibration. The results of the verification step are compared with an observed data set, specifically 
reserved for this purpose. 
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In case of hydraulic and hydrological models, where time is involved, model calibration is usually done 
in three steps: 


(1) Selection of a simulation time period from the available observation data set. Normally the 
available data set is split in two, one set for calibration and one for verification. Selection of the 
time length for calibration depends on the problem that needs to be addressed by the model. For 
example if flood is of interest the observed period should contain a flood event. 

(2) Preliminary calibration, when based on the experience of the modeler a set of parameters are 
chosen and changed in such a way that the simulation results are as close as possible to the 
observed data. 

(3) Refined calibration, when a thorough analysis of the output results is done and refinement of 
parameters values is carried out (either manually or by automatic procedures). 


In order to determine the quality of the calibration a set of criteria are chosen, based on which the fit 
between simulated and observed data is evaluated. The selection of criteria depends on the objective of 
the modelling. In flood event modelling, for example several criteria such as time to peak, peak of the 
hydrograph, are used to assess the correspondence between observed and simulated flows. 

Parameter calibration is not an easy task because some parameters may be influencing each other 
(especially in complex models) and this can generate sub-optimum determination of parameters. It is 
advisable that some of the parameters are obtained by direct field measurement, where possible, and 
for models for which they have physical meaning. Bathurst (1986) advises that a sensitivity of the fitting 
criterion to any changes in parameter values should be carried out before selecting any value for a 
parameter. 

Immediately after calibration verification is carried out, by using the remaining data sets from the 
available observations. While verifying the outputs of a simulation the assumption is made that the 
governing equations and computer coding are valid, because the model passed the validation test. 
Verification therefore needs to be objective, subject to rigorous testing criteria, while taking into account 
the conditions for testing (e.g., years, inflow flood hydrograph, land use, etc). 

With the development and application of models a special issue arises, the uncertainty of the model 
results. Uncertainty in modelling is of three types; structural, parametric or stochastic (Smith & Kuhn, 
2002). Structural uncertainty comes from lack of knowledge about the phenomena that is modeled (i.e., 
lack of measurements, poor measured data, human errors, etc), while parametric uncertainty comes from 
using parameters with unsure values. Stochastic uncertainty shows the randomness of the phenomena and 
can not be reduced. 

Though uncertainties are important in modelling, models can be evaluated based on clear defined 
criteria. 


1.4 DEVELOPMENT AND APPLICATION OF MODELS 


As a summary of the notions presented in this chapter a framework for development and use of models in 
engineering practice is presented bellow and schematized in Figure 1.7. 

Several authors (Branson et al., 1981; Shepherd & Geter, 1995) have outlined the steps that need to be 
taken in order to develop a proper modelling approach for solving waer related problems. These steps are: 


(1) Define the objectives of the model, which requires that the modeler addresses the issues that are 
posed by the problem to be solved and determines what is the type of problem (water quantity, 
water quality, flood volume, etc); 
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Figure 1.7 Model development and application steps. 


Q) 


(3) 
(4) 
(6) 
(6) 


(7) 
(8) 
(9) 
(10) 


Determine if there is available a similar model that can be used. If there is no model available than 
anew model is defined and implemented on a computer; 

Define the mathematical representation of the phenomena that requires solution; 

Select the numerical approach that solves the mathematical representation of the phenomena; 
Coding and implementing the code on a computer; 

Validation of the mathematical model, by checking the model against simple examples. There is 
a need to check what are the limitations in applicability, and what are the required boundary and 
initial conditions; 

Sensitivity Analysis for the parameters of the model; 

Testing and Evaluation, when calibration and verification is performed; 

Application of the model for different problems; 

Presentation of results in graphical format, tabulated format or animated form to decision makers 
and stakeholders. 


A good modelling practice does not mean just the development of a model and its implementation on a 
computer, but also continuous maintenance and refinement (improvement) without affecting its integrity. 
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Chapter 2 
Modelling water related problems 


2.1 BASIC CONSERVATION EQUATIONS 


Computational fluid dynamics (CFD) is an important tool in engineering problem solving practice, covering 
a whole range of disciplines. As described in the previous chapter, Hydraulics is just one of the application 
fields of CFD, and it encompasses the study of free-surface and pressurized water flows. Traditionally 
Hydraulics is the scientific basis for engineering applications (Jain, 2001; Rousse, 1950). Hydraulics refer 
to water in both its states; still and flowing. The branch of Hydraulics concerned with the study of water 
flowing is referred as Hydrodynamics or Fluid Dynamics, when the studied fluid has different properties 
than water. The scope of the book is to address the equations representing the flow of water. 

Conceptual models, based on empirical relations tested through observations (field measurements or 
laboratory experiments) are the first types of mathematical representations that were used in Hydraulics. 
Empirical relations were not the only representations used to describe behaviour of a water element; 
mathematical formulations in form of equations, were also used, however these equations could only 
give descriptions of the water behaviour and were representing relations between different variables 
and properties of the water. These types of equations were too complex to be solved analytically. The 
development of computers allowed for the formulation of solutions of particular situations, which brought 
the understanding of essential features of phenomena. 

In this chapter the general formulation of the equations of water related problems are formulated. Their 
solution approach using numerical approximation is detailed in the following chapters of the book. The 
general classification of the equations, from mathematical point of view, is also given, because of the 
specific properties that each of the equations has. 

In continuum mechanics the equations representing the general behaviour of a continuum are referred 
to as conservation laws, and they are based on Newton's laws. In Hydraulics the laws are conservation of 
mass, energy and momentum equation. These laws are not applicable only to fluids, they are also valid 
for solids, however they differ in behaviour, therefore due to the scope of the book, only the water related 
equations are addressed. 

The laws of physics that are used for determining the general principle of conservation laws are 
Newton's laws of motion and the first and second laws of thermodynamics. 
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Newton's laws can be stated as: 


* First law: a body is at rest or continues in a uniform state of motion, unless acted on by an applied 
force. 

* Second law: the net force acting on a given mass is proportional to the time rate of change of the 
product of mass and velocity (also known as linear momentum); 

* Third law: action and reaction are equal and opposite. 


It is the second Newton’s law that introduces the notion of change of state. 


The laws of thermodynamics can be stated as: 


* First law: in a system thermally isolated by impermeable walls the work done in taking one state A 
of a system into another state B is entirely determined by the terminal states A and B. The internal 
energy difference between state A and state B is defined as the mechanical work done in taking 
either state A into state B or state B into state A; 

* Second law: there is a tendency on the part of nature to proceed towards a state of greater disorder. 


Based on the above laws the general principle of conservation laws is that the rate of change of a 
quantity u within a volume V plus the flux of u through the boundary A (noted as f(u)) is the same as the 
rate of production of u denoted by S(u, f) (see Figure 2.1): 


2 usav s f fonda- f stunav =0 (2.1) 


fa) 


Figure 2.1 Control volume of quantity u in a (x, y, z) referential system. 


Equation (2.1) is referred to as the integral form of the conservation law, and can be further detailed for 
mass, energy and momentum. 


2.1.1 Conservation of mass 


Conservation of mass states that for any control volume, during a small time interval Af the mass 
entering the volume minus the mass leaving the volume equals the change of mass inside the control 
volume. In case of continuation of mass equation (2.1) for mass m of density p, and advection 
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velocity u over the control volume V gives: 


Sf p-av+f oin dn co Q2) 
where: 
fu)=p-u and u-n=u (2.3) 


In case of conservation of mass the term S, is the source term, and since there is no mass production, S 
is zero. Integration of equation (2.2) yields: 


Ó 8 


EC NCNCA RR 
with V — ( mm 2) the so-called nabla operator. 


Equation (2.4) can be further detailed if the product rule for vectors is applied. It yields: 
Pg 7 2.5 
3: tE P + pV-u=0 Q.5) 


which can be re-arranged as: 


Dp = 
ix 2.6 
Di + pV -u=0 (2.6) 


where D/Dt is the total derivative with respect to time: 


D ə dð dyd doa 
Dt ot dtox dt dy dt az 


In general water is considered an incompressible fluid (i.e., p = const.), hence the continuity equation 
simplifies to: 


V-a=0 (2.8) 


Equation (2.6) is the continuity of mass, while equation (2.8) can be regarded as the equation of 
continuity of volume. 


2.1.2 Conservation of momentum 


Conservation of momentum is demonstrated in x direction only. Conservation in y and z dimensions is 
similar. 

In case of momentum in x direction, the conserved quantity is pu,. In equation (2.1) momentum flux in 
x direction is pu,u and S are the body and boundary forces acting on the control volume. Body force is the 
gravity and boundary forces are pressure, shear and surface forces. 
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Momentum equation in x direction, according to equation (2.1), is: 


o(pu,) " 
UV - i 
gp a NUM. Q3) 
Pressure force | Friction Force Gravity Force 
If gravity and pressure forces are detailed, yields: 
D op 
udo 2.10 
D 3x TPE t+ Sra Q.10) 
In three dimensions equation (2.10) becomes: 
Du = 
—=-Vp+ pg + 2.11 
P p+ pg * $, (2.11) 
where i = (u,,u,,u,) and Ẹ = (8,.8,.8.)- 
Using the total derivative definition, equation (2.11) becomes: 
Apri) E E EE EA (2.12) 
In equation (2.12) the term u x u is the cross product, introduced, in calculus, by definition, as: 
Ha. UU, uu 
uxu=Vi uu, UU, UU, (2.13) 
UU, UU, Uu 
Hence, 
pu,u, puu, puu, 
V(pu x ü) = V| puu, pu,u, puu, 
pu.u, puu, pu.u. 
2 (pu) + S ta) + 2 pua) 
(2.14) 


2 (pui,) + S ouu) + 2 upu) 


d Ó d 
| ax Ms) T ay Pd) + 3; Pete) | 


- (Vp-i)u + pu(V -u)+ pu: V)u 
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Left side of equation (2.12) becomes: 


MOD opis) 2 pos iP sawp. ii) + pi(V-ü) + pli- Vi 


ot (2.15) 


E op a . . . 
(pev gio |-0-continuity equation 


Based on (2.15) equation (2.12) becomes: 


à L . 8 - 
PA + pū- V) = -Vp- pg + S, (2.16) 
Equation (2.16) represents the conservation of momentum in a three dimensional (3D) space. 


2.1.3 Conservation of energy 


For phenomena in which temperature varies, in addition to momentum equation, conservation of energy 
must be considered. The energy equation accounts for both kinetic and potential energy. Hence energy is 


expressed by the term ple + i. and the equation for conservation of energy is: 
d u? w ). Š = i 
g Pl et a +V. petu = V.(KVT) - V - (up) + V - (ut) + pgu (2.17) 


where e is the potential energy; p is pressure; 7 is stress tensor; k coefficient; and T is temperature. 

The three presented conservation laws form the so-called Navier-Stokes equations, which are space 
and time-dependent, that is, there are four independent variables. In a three dimensional domain there are 
six dependent variables and only a set of five equations available; a continuity equation for conservation of 
mass, three equations for conservation of momentum and one equation for conservation of energy. The six 
dependent variables are pressure, density, temperature, and three components of the velocity vector u. The 
six equation that is used is the law of thermodynamics relating pressure and temperature: 


=p: V = nRT (2.18) 


where p-pressure; V-control volume; T temperature; and n, R coefficients. 

A system of six equations with six unknowns is formed, which can be solved. Usually, the Navier- 
Stokes equations are too complicated to be solved, and simpler forms are used. These are discussed further 
in this chapter. 


2.2 MATHEMATICAL CLASSIFICATION OF FLOW EQUATIONS 


One of the conclusions of the previous paragraph is that the physical processes related to water problems 
are represented mathematically by differential equations (commonly abbreviated to ODEs or PDEs, for 
Ordinary Differential Equations and Partial Differential Equations, respectively). A differential equation 
is a relationship among various distinct derivatives of a function of one or more variables. Examples of 
such processes are propagation of waves in a river channel, pollutant transport, reservoir operations, 
and so on. In describing these phenomena rate of change of a variable is evaluated, which results in the 
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representation of the physical processes in terms of differential equations (DE). Differential equations 
are models of the phenomena and require integration to be solved. This is not always possible, hence they 
are solved numerically. 

As opposed to an algebraic equation, where the unknown is a number, in a differential equation, the 
unknown is a function that expresses the behaviour of the variable as a function of time and space. A 
differential equation has order, which is given by the order of the derivative of the unknown function 
involved in the equation. If the values of the unknown function and its derivatives at some point are known 
the solution to DE is called an initial value problem (IVP). If no initial conditions are given, the description 
of all solutions to the DE is called the general solution. 

DE are classified taking into account several criteria, as follows: 


* dimension of the unknown function: 
o ordinary differential equation (ODE), when the unknown is a function of one independent 
variable only, for example, u(t); 
o partial differential equation (PDE), when the unknown is a function of multiple independent 
variables, for example, u(t, x, y); 
* number of equations: 
o single differential equation; or 
o system of differential equations (coupled); 
* order: 
o the highest order of the derivative that appears into equation, that is, an equation is of n-order if 
the DE has a n-th derivative, and not higher; 
* linearity: 
o linear, when all terms are linear in unknown function and its derivatives; 
o non-linear, when the linearity does not hold. 


In water related problems, most frequently the 2nd order, linear PDEs with two independent variables 
(homogeneous or nonhomogeneous) are used. In such a case, if the two independent variables are assumed 
to be x and f, the PDE is written as: 


ou 
or? 


= Dan + Eon% + F(x,t)u(x,t) + G(x,t) (2.19) 


o? o? 
dp s y BG. t Cf) 


where the coefficients A(x, f)... G(x, f) may be functions of x, t, or both, or they may be a constant. The 
function u(x, f) is the solution of the equation. The second-order PDE (2.19) has a discriminant defined as: 


A = [BG OP - 4: AGsf) : CCo f) (2.20) 


Depending on the sign of the discriminant A (as defined by equation (2.20)), equation (2.19) is 
classified as: 


* elliptic PDEs, if A < 0, which are PDEs with smooth solutions, easy to solve; 

* hyperbolic PDEs, if A » 0, which have solutions that may have discontinuities, usually difficult to 
solve computationally; and 

e parabolic PDEs, if A = 0, which may have features of hyperbolic and elliptic PDEs combined. 
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The terminology for this classification is coined from the geometry of conic sections that satisfies a 
similar type of equation like (2.19), just that they are expressed in coordinates x and y. 

Each type of equation describes a certain type of physical phenomena; elliptic PDEs describe processes 
that are in steady state, time-independent (e.g., steady-state aquifer flow, backwater curves); hyperbolic 
PDEs describe time-dependent, conservative physical processes, (e.g., advection/convection), that are 
in unsteady state and not evolving towards steady state; and parabolic PDEs describe time-dependent, 
dissipative physical processes (e.g., such as diffusion, transient aquifer flow) that are evolving towards a 
steady state. These three types of equations are very important for flow problems and the solution for each 
one of them is further presented. 


2.2.4 Solutions of ODE 


An Ordinary Differential Equation (ODE) is a differential equation where the solution depends on one 
independent variable only. The general form of an ODE is: 


du(t) _ 2.21 
7 7 ft) up 


where u(f) is an unknown function of the variable f (time) and fis an arbitrary function in u(f) and t. 
Ordinary differential equations may have analytical solutions that are easy to obtain if the function 
flu, t) is not too complicated. However even if the analytical solution exist the implementation of the 
solution into a computer programme may be difficult. 
If along with equation (2.21) an initial value of the function u(x, f) is given then the problem is called 
initial value problem (IVP) and can be solved easily using numerical approximations, as it is detailed 
further in Chapters 4 and 5 of the book. 


2.2.2 Solutions of PDE 


In order to have a solution for a PDE, the problem must be well posed mathematically (Garabedian, 
1966; Cunge & Verwey, 1986; Ligget & Cunge, 1975). This implies that solution fulfils three conditions; 
it exists; it is unique; and depends on auxiliary conditions such as boundary conditions. Fletcher (1998) 
states that uniqueness of solution is not a problem in general, and if it does not exist is mainly due to 
failure to fulfil auxiliary conditions. Moreover from mathematical point of view for some particular 
equations there are isolated points in the definition domain for which the solution may not be unique. 
There are physical phenomena for which multiple solutions exists, due to the physical phenomena, not 
due to the mathematical position of the problem. Such an example is the case of transition from laminar 
to turbulent flow. 

Any PDE can be solved if boundary conditions and initial conditions on the computational domain are 
known. Typical representations of the computational domains and required boundary conditions for the 
three defined types of second order PDEs are represented in Figure 2.2. 

The definition of boundary conditions constitutes the start for determining the solution of any problem 
inside the computational domain. There are three ways of defining boundary conditions; Dirichlet type 
of conditions, when the values of the unknown function u are known at the border of the computational 
volume (on A, see Figure 2.1); Neumann type of conditions, when derivatives of the unknown function 
is known at the border of the computational domain; and Robin or mixed type of conditions, when a 
combination of Dirichlet and Neumann conditions are applied at the domain boundary. 
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Figure 2.2 Computational domains, boundary and initial conditions for PDEs. (B.C. Boundary conditions; 
I.C. Initial conditions; (a) Hyperbolic type of equations; (b) Parabolic type of equations; (c) Elliptic type of 
equations). 


An important issue in finding the solution of a PDE are the errors. Hyperbolic types of problems are the 
most exposed to errors because in their case conditions at the boundary of the computational domain are 
the ones introducing the errors that will propagate inside the computational domain. 


2.2.2.1 Boundary conditions for fluid flows PDEs 


Majority of the flow phenomena in hydraulics (hydrodynamics) are governed by a system of two 
differential equations (as it is for example the case of free surface flows). A system of two first order 
partial differential equations, having two independent variables (€, 7) and two unknown functions 
(u and v), is written as: 


Ou ou ov ov 

atata hon = f (2.22a) 
Ou ou ov ov 

a seth a tose tha = fh (2.22b) 


where, ó and 77 are the independent variables. For water related problems ó corresponds to the x-coordinate 
and 7] to time t. The unknown functions are variables that correspond to flow quantities, such as flow 
velocity or water depth h. 
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Because both u and v are functions of (E, 1), their total derivatives in the (€, 7) plane are: 


du = aE F an (2.23a) 
ov ov 
dy = prés t 25 4 (2.23b) 


These derivatives can have different values in different regions of the (é, 7) plane, or may not exist in 
some regions. The curves that are splitting the lines into different regions are called characteristic curves 
and are determined from the condition that the determinant of the system formed from the four equations 
(2.22a, b) and (2.2.3a, b) is not zero, that is, the system of equations has a unique solution. The formed 
system of equations is written in matrix format as: 


a b c d || duldg fi 
a b, c, d,|Oulm| |f (2.24) 
dE dn 0 Oj|owoé| | du 
0 0 dé dm|owom| |dv 


The unknowns of the system (2.24) are the derivatives of the unknown functions u and v. These are not 
determined if the determinant of the system (2.24) is zero, The value zero of the determinant defined by 
(2.24) determines the regions in space beyond which there is no solution for the given equation, because 
the derivatives do not exists. The condition is written as: 


a b c d, 
b d 
- o a (2.25) 
dé dn 0 0 
0 0 dé dm 


Calculating the determinant gives: 
dn Y dn 
(ac, — aci) - (2 — (ad; — ad, + bc, — be): dé + (bd, — bad) = 0 (2.26) 


Rearranging yields: 


any a, ,. 
«(22 TE ere (2.27) 


with coefficients: 


a = (ac, —a,¢,); b-(-ajd,*ajd — bic + bac), c= (bid, — badi) (2.28) 
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Equation (2.27) has two solutions, real or imaginary, depending on the value of the discriminant. The 
two solutions of the equation (2.27) are the known solutions of the quadratic equation: 


bua» - 
dn _ —bt«b^-4ac (2.29) 


dë 2a 


Equation (2.29) represents curves in the (5, 7) plane, called characteristic curves (C1 and C2), and 
determines the number of required boundary conditions in order to have a unique solution. The discriminant 
is similar to the one defining the types of PDEs, hence depending on the value of the discriminant the 
problem can be of 3 types: hyperbolic — with two real characteristics; parabolic — with one characteristic; 
or elliptic — with two imaginary characteristics (see Figure 2.3). 


(a) ; Zone of (b) " 


influence (c) t 


C1-2C2 


Domain of 
dependance 


Zone of 
influence 


Hyperbolic P.D.E. Parabolic P.D.E. Elliptic P.D.E. 
Figure 2.3 Domain of solution for a system of two first order differential equations. 


If the analysis is done from the point of view of the solution of flow phenomena it can be stated, as 
already mentioned, that hyperbolic behaviour means that a system develop from the state at point P of the 
computational domain towards a state in between the characteristics C1 and C2 (shock waves). In case of 
parabolic type of problems it can be seen that from a state P in the computational domain the solution can 
evolve only in one part of the computational domain. Elliptic type of problems can develop from any state 
P to any state in the computational domain. 


2.2.2.2 Hyperbolic equation examples 


The typical example of a hyperbolic equation is the wave equation: 


Qu  ,0O'u 


DU gee scs 2.30 
an ^ a 2-30) 
where u is the unknown function and a is the wave speed, considered a constant in this example. 
Given the initial and boundary conditions, over the domain (0, L) (see Figure 2.42): 
u(x,0) = sinzx 
ou =0 (2.31) 
(x-0) 


u(0,t) = u(L,t)= 0 
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an exact solution can be determined as: 
u(x,t) = a: sinZx- coszt (2.32) 


Representation in the (x, £) plane, of the solution given by (2.32), is shown in Figure 2.4b. 


(a) Initial Condition (b) 


Sin(x*x) 
Solution 


0 0.2 0.4 0.6 0.8 1 


Figure 2.4 (a) Initial condition and (b) solution of the wave equation. 


In general, for different types of initial conditions the solution of the equation (2.30) is not easy to 
be determined. The method used to solve the wave equation is usually the method of characteristics 
(Abbott & Basco, 1982). According to equation (2.27) there are two real characteristics available for 
equation (2.30): 


dt 

2.33 
dsl (2.33) 
dt 


The set of equations (2.33) are two straight lines (because a is a constant), which show that a disturbance 
in solution of u, taking place at point P influences the solution domain after P (CPD in Figure 2.5(b)). 
Solution in point P is influenced by disturbances in the domain before point P (APB). Depending on the 
position of P in the computational domain initial conditions may be sufficient to determine the solution 
at point P, uniquely. However there are points in the computational domain for which solution can be 
determined uniquely only if boundary conditions are given (see Figure 2.5). 

Another example of hyperbolic equation, which is very common for many phenomena in aquatic 
environment, is the advection equation: 


du ao (2.34) 
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In case of equation (2.34), just one characteristic line passes through every point of the space (see 
Figure 2.6). Equations (2.30) and (2.34), presented herein, are equation showing propagation phenomena. 
No dissipation effect is present, just convection. 


(a) C Characteristics C* Characteristics (b) 
[n ^ [a Tw t D C 
Domain of 
influence 
P / 
f, / 
YY YY mae of 
ependance 
A B C - 
Specified conditions A B x 
X Wave equation 


Hyperbolic equations 


Figure 2.5 Characteristic curves for hyperbolic PDEs. 


ta 


Figure 2.6 Characteristic lines for the advection equation. 


2.2.2.3 Parabolic equation example 


Phenomena in nature for which dissipation mechanism appears along with propagation are of parabolic 
type. The diffusion of a pollutant in a water body is such an example. The diffusion equation in one 
dimension of space is: 


ou e pou (2.35) 
X 


where D is the diffusion coefficient. Equation (2.35) is defined over the domain (0, L) (see Figure 2.7). 
Given the same initial and boundary conditions as in the case of wave equation (2.30), the analytical 
solution is: 


u(x,t) = sin mx - exp?t) (2.36) 


The exponential part of the solution given by equation (2.36) represents the time decay. 
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Initial Condition 


Sin(n*x) 
Solution 


Figure 2.7 (a) Initial conditions and (b) solution for the diffusion equation. 


When other types of boundary and initial conditions are given the solution has to be computed either 
through integration or through numerical approximation. In case of parabolic equations characteristics are 
not so important because they are equal, and they show that solution at point P of the computational space 
can influence any point in the computational domain. In case of parabolic equations any combination of 
Dirichlet, Neumann or Robin boundary conditions are appropriate to be used. 


2.2.2.4 Elliptic equation example 


Elliptic equations are time independent and therefore only valid in case of steady flows. A typical example 
of elliptic equation is the Laplace equation: 


Qu Qu 

L^ 44>" =90 (2.37) 

ox? dy? 

In case of elliptic PDEs all boundary conditions should be given (see Figure 2.8). The characteristic 
lines are of complex nature and cannot be represented in the real domain, hence they are not used. The 
important characteristic of an elliptic equation is that any disturbance at point P in the computational 
domain will influence all the points in the computational domain. 


Computational 
Domain 


B.C. 


B.C. * 


Figure 2.8 Computational domain and boundary conditions for elliptic PDEs. 
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2.3 NAVIER-STOKES AND SAINT-VENANT EQUATIONS 
2.3.1 Navier-Stokes equations 


The dynamics of fluid flow is governed by the Navier-Stokes equations, which are derived from the 
conservation laws (mass, momentum and energy), as shown in Section 2.1 of this chapter. Derivation of 
Navier-Stokes is shown in detail in many text books (Hirsch, 1988; Popescu, 2013). In case that flow of 
water is studied, several assumptions are made on the fluid; it is incompressible; viscous; homogenous; and 
does not conduct heat. The obtained equations are called incompressible Navier-Stokes equations and are 
written as a system of equations of continuity and momentum (Popescu, 2013). In case of incompressible 
flows equations (2.8) and (2.15) are: 


MCCC: E 
gy = Qc ur Vu FV +S (2.38) 


Vu =0 


where u is the velocity vector, Ù viscosity (controls the velocity of diffusion), p is pressure and S the 
external forces. First equation of the system (2.38) shows the change in velocity, which is composed of 
a summation of four terms. The first term in the right-hand side of the first equation of system (2.38) is 
called advection (convective) term, and it shows how the fluid moves around, it shows the force exerted on 
a particle of the fluid by the other particles of fluid surrounding it. The second term in the same equation 
is the diffusion term and it shows how the fluid motion is damped. The third term is the pressure term, and 
the fourth term are the external forcing functions such as gravity or wind. 

Navier-Stokes equations are not easy to solve, because of the non-linearity of the advective term of the 
equation. Analytical solutions however, are available for particular cases, very simple ones that have a 
simple geometrical computational domain and simple boundary conditions. 

For flow of water the Navier-Stokes equations are solved, by making further assumptions, which simplify 
the equations and allows for numerical approximation of the solution. Each of the simplified approach is 
used for different problems, because the assumptions made are depending on the considered phenomena. 
There are three main approaches, in water related flows, in order to represent Navier-Stokes equations in a 
simpler form: 


e Averaging the equations: Reynolds averaged equations plus a turbulence model (RANS); 

* Large Eddy Simulation (LES); 

e Approximations of the averaged equations, that is, equations averaged over the cross-section (Depth 
averaged). 


In the first approach, the average of the flow is used, and the turbulence fluctuations are considered less 
important. Each physical quantity is represented as an average over the computational time interval. The 
quantity itself is represented as a summation between a mean value and a fluctuating part. The averaging 
transforms the non-linear terms in linear ones, but they do introduce a new term, the derivative of the 
averaged quantity. The equations obtained are called Reynolds averaged equations. The additional term 
in the equation requires the use of a turbulence model. The LES forms of the Navier-Stokes equations 
are used for turbulence modelling (Aldama, 1990; Lesieur, 1997; Sagaut, 2001). In the third approach 
depth averaged approximations of velocity in space are used for the simplification of the Navier-Stokes 
equations. 
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2.3.2 Saint-Venant equations 


The flow of water can be in a closed or open space. These two types of flow differ in the fact that in the 
open space water flows with a free surface, while in the closed space is a pressurised flow. The flow in a 
closed space is flow in conduits (pipes) and can be in some instances with free surface, however they can 
be considered the same as the flow in an open space. The free surface flow is the flow occurring in oceans, 
estuaries, coastal regions, lakes and impoundments. This kind of flow is subject to atmospheric pressure 
and water surface changes in space and time, hence it is more difficult to solve. Variables that are part of 
the equation; depth of flow, velocity or discharge, bottom slope, position of the free surface flow; are all 
interrelated. Moreover in case that flow is in a river the shape of the cross-section of the channel flow is 
arbitrary, while in pipe flow cross-sections are simple (circular or regular shape). Roughness of a channel 
may vary much more than the roughness of a pipe, which makes the free surface flow more complex. 
These equations are not easy to solve, hence research on methods of solution of the free surface flow 
equations has received considerable attention. 

Due the complexity and importance of free surface flow, the equations are presented and analysed 
in detail here. Methods of solution for these equations through numerical approximations are shown in 
Chapters 4 and 5. For the special case of flow in rivers and lakes these equations and their solution is 
addressed in Chapters 7 and 8 of the book. 

Free surface flow is unsteady in nature and the equations expressing it are unsteady flow equations. 
There are however cases in which the problems are steady in nature and the solution is easier to be 
determined. Though the equations of steady flow are determined separately, through assumptions on the 
steady character of the flow, they can also be determined from the unsteady flow equations, as particular 
case, when all changes with time are set to zero. 

Flow is classified in respect with its variation with space, not only with time, as uniform and non- 
uniform flow. Flow is considered uniform if depth and velocity of flow is the same at every point in space. 
It occurs when gravity forces are in equilibrium with friction forces. 

In case that the wave length of the free surface flow is much more bigger than the depth of flow then 
the equations describing the phenomena are referred as shallow water flow equations or Saint-Venant 
equations. These equations are Reynolds depth averaged forms of the Navier-Stokes equations. They can 
be determined also independently of Navier-Stokes equations, using physical arguments and by making 
assumptions on the fluid flow. 

For simplicity of formulation Saint-Venant equations are further discussed for the case of one dimensional 
flow. The basic assumptions for the derivation of Saint-Venant equations are (Saint-Venant, 1871): 


(1) The wave length of the disturbance of the flow is very long as compared with the depth of the flow. 
This assumption implies that flow is one dimensional and parallel to the walls and bottom forming 
a channel in which flows takes place; 

(2) The channel geometry is not changing in time, that is, deposition or scour of sediment is very small 
and considered to have no effect; 

(3) Streamline curvature of the channel is small; lateral and vertical accelerations are negligible 
compared to the longitudinal accelerations; therefore, the pressure distribution is hydrostatic; 

(4) Slope of the channel bed is small, that is, the tangent and sine of the angle between bottom and 
horizontal have the same value as the angle itself; 

(5) The effect of boundary friction force is estimated using a relation derived from steady uniform 
flow; 

(6) Water surface in any cross section of the channel is assumed to be horizontal and parallel with the 
bottom of the channel. 
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There are two ways of expressing the Saint-Venant equations; in conservative form, when the dependent 
variables are discharge and depth of flow; and in non-conservative form, when velocities and water depth 
are the dependent variables. Equations expressed in velocities are non-conservative because velocities are 
not subject to conservation, and the equations do not hold across discontinuities, such as hydraulic jumps 
or in case of shocks. Equations in conservative form have the advantage that they do hold even if some of 
the assumptions of Saint-Venant equations break, as it is the case for example of the hydraulic jump. They 
are used when there are significant variations of flow in space. 

The Saint-Venant equations are obtained by writing the conservation of mass and momentum over a 
small control volume. For a 1D channel, which has the longitudinal profile, and cross-section as the one 
represented in Figure 2.9, Saint-Venant equations for an arbitrary cross-section of a channel, in conser- 
vative form, are (Cunge et al., 1986; Popescu, 2013): 


0A . a0 
= 2; 
d ax ^ use 


2 
dQ , 20 00 +( Q rs JA a dS, cs d (2.39b) 


ot A ox A> B ) ox 


md | Kinematic wave 
——————————————ÁÁ—— | Diffusion wave 
————————————————————À—À———— | Dynamic wave 


h 

Figure 2.9 Control volume of a continuum fluid element. 
Saint-Venant equations in non-conservative form in 1D are: 
ðh | O(uh) 
CEA e - 2.39c 
ot ox ( ) 
d £u g 9$, + gAS, = qu (2.39d) 
ot ox d f f 


————————————— | Kinematic wave 
— ———————————8———— — ——— —— | Diffusion wave 
eee tea eli ee. eee tou Mee Dni ete uc. ian eicere | Dynamic wave 
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In equations (2.39a—d) the notations are, as per Figure 2.9; A, area of the cross-section; Q discharge in 
the control volume; q lateral flow; A water depth; S, the slope of the channel, S, the friction slope; u velocity 
of flow; and u, velocity of the lateral flow. 

Equation (2.39a) and (2.39c) represent the continuity equation, while (2.395) and (2.394) represent the 
momentum equation. 

In both representations, momentum equation has five terms; the local acceleration, which describes 
the change of momentum due to the change of velocity over time; the convective acceleration that shows 
the change in momentum due to change of velocity along the channel; the pressure term that describes 
the pressure force which is proportional to the change in water depth along the channel; the gravity term 
that shows the force due to gravity; and the friction term that shows the friction force proportional to the 
friction slope. 

The system formed by the above two equation, in conservative or non-conservative form, does not have 
an analytical solution. There are cases when momentum equation can be used with less than five terms, 
depending on the problem to be solved. Some of the terms, of the momentum equation might be very 
small, as compared to the other terms, hence they can be neglected (Price, 1985). Depending on how many 
terms are used in the representation of the equation of momentum, the following terminology is used: 


e Kinematic wave representation, when momentum equation reduces to two terms: S} = S; 

* Diffusive wave approximation, when the first two terms of the momentum equation, called inertia 
terms, are neglected; and 

* Fully dynamic wave, when the full momentum equation is considered. 


Solution approaches for these three cases are detailed below. 


2.3.2.1 Kinematic wave solution 


In the case of a kinematic wave approximation, the acceleration and pressure terms in the momentum 
equation are neglected. Kinematic wave approximation has no terms varying in time, the terms in the 
momentum equation represent the steady uniform flow. However the flow itself is not steady, because the 
unsteadiness is taken into account by the continuity equation. In case that there is no lateral flow (q, = 0), 
the system of equations (2.39) has an analytical solution. 

The slope of the energy line is the same as the slope of the channel bed, which is equivalent with stating 
that discharge is a function of flow depth, Q = Q(ARG9), hence: 


dQ _ dQ dh 
m ae ae (2.40) 
Inserting (2.40) into (2.38) it yields: 

dA dQ 90 _ 9A 90, AQ dh _ on 


dO g^ oc 00 of Uh ax | 


Equation (2.41) is a form of advection equation, with the wave celerity c = 2o Re-arranging (2.41), for 
the case of zero lateral flow, yields: 


9A dQ dh, dQ Oh - 
dQ dh Ə dh Ox 


0 (2.42) 
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Hence: 

oh oh 

ee PELLE. 2.43 
a To X 0 (2.43) 


An expression for the kinematic wave celerity c may be obtained from one of the steady flow equations 
such as Chezy or Manning (Chow er al., 1988). For example by using Manning equation, and a wide cross- 
section approximation, for which hydraulic radius is approximately the same as the water depth, the value 
of the wave celerity is: 


fF (2.44) 


Equation (2.43) is the kinematic wave equation and is a linear advection equation if c is a constant. The 
advection equation with a positive velocity c will propagate just in the positive direction of the x axis, that 
is, towards downstream. In general, the wave velocity is a function of water depth c = c(h). The bigger the 
flow depth, the bigger the value for c. For one-dimensional problems this results in a steepening of the 
wave front. Kinematic wave equation applies to situations where flows are on steep bed channels, but do 
not apply to backwater effects, nor to tidal flows. An extensive study on the kinematic wave can be found 
in Ponce and Simon (1977). 


2.3.2.2 Diffusive wave solution 


The kinematic wave model is suitable just for some problems as it is for example the case of translation of 
a flood wave. A second approximation of the momentum equation is the one in which the pressure term is 
taken into consideration as well, hence: 


OF oe E. (2.45) 


which is equivalent with stating that O = Of n(x), 94), Replacing this relation into equation (2.38a), and 
considering lateral flow to be zero, lead to: 


2 
g-o a 


where c is the kinematic wave celerity (as previously defined) and D is the diffusion coefficient. Equation 
(2.46) has the form of an advection-diffusion equation, hence its name. The diffusion coefficient has the 
expression: 


X em 
. B |.( 0h 

(as) 
The diffusion coefficient is always positive because the energy slope is always positive. Different from 


kinematic wave, in case of diffusive wave, as the name is stating, there is diffusion. While the water wave 
travels towards the downstream direction diffusion in form of attenuation happens. 


(2.47) 
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Applying the diffusive wave description of phenomena it is possible to solve problems with relatively 
small backwater effects, slowly propagating flood waves and runoff on mild slopes. Several methods of 
solution are available in the literature for the diffusive wave approximation such as Muskingum-Cunge 
method (Cunge, 1969), the Zero-Inertial solution (Strelkoff & Katopodes, 1977) or Parabolic and Backwater 
(Todini & Bossi, 1986). 


2.3.2.3 Fully dynamic wave solution 


In the case that inertial forces are as important as the pressure forces the fully dynamic expression of the 
momentum equation is used for the solution of Saint-Venant equations. Due to their complex nature these 
equations can not be solved analytically, however they can be solved using numerical methods, that is, 
using computer power to solve a system of equation obtained after the application of different numerical 
schemes to equations (2.38) or (2.39). 

Popescu (2013) makes an overview of the different numerical techniques available for solving the Saint- 
Venant equations. These methods include the semi-analytical solution by the method of characteristics 
(Abbot & Basco, 1998), explicit difference methods, semi-implicit methods (Casulli, 1990), fully implicit 
methods, and Godunov methods (van Leer, 1979). 

Solving Saint-Venant equations numerically implies that conditions on stability and convergence of 
the numerical solution must be met. Errors are mostly introduced by explicit methods and are avoided 
by imposing the Courant Frederichs Levy (C.F.L.) condition (At € Ax/u). Implicit methods do not 
induce stability problems, but are subject to constraints in the selection of the time and space step of 
computation. 


2.3.3 Characteristic form of Saint-Venant equations 

Saint Venant equations, as they are formulated in (2.39), can be transformed into the so-called 
characteristic form, which is useful to understand their meaning. Along with the formulation (2.39), the 
auxiliary variable c = ,/gh is introduced, for which the derivatives with respect to independent variables 


x and f are: 
oh 2c dc 
oh 2c Oc 
27 d pm (2.48b) 


Equations (2.48) are substituted in (2.39a) and (2.39b), then the obtained equations are divided by c 
and multiplied by g. As a result, the following two equations are obtained: 


e ay ae eee) 
ou ou dc 
3r +u- zte J =0 (2.49D) 
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Adding (2.49a) and (2.49D) and then subtracting them, results in: 


HUF 29) , qu 0) 90120 ig (2.504) 
(u — 2c) O(u — 2c) _ 
= bipes ~~ = 0 (2.50D) 


Recall that the total derivative of an arbitrary function f(x, f) is: 


df of dx of 
- : 2.51 
dt ot dt ox ee 
dx _ df _ C : . ] 
If atte then di 0. This is similar with saying that: 
d(u + 2c) dx 
= E e 2.52 
di 0 along di (u+c) ( a) 
or equivalent 
dx 
J, = (u+ 2c) 2 const along a^ (u *- c) (2.52b) 
Similarly: 
d(u — 2c) . dx 
= D ees 2.53 
di O along line di (u — c) ( a) 
or equivalent 
dx 
J_ = (u— 2c) = const along m (u — c) (2.53b) 


Equations (2.52a,b) and (2.53a,b) are the characteristic equations; J , J, are called Riemann invariants. 
These equations are of hyperbolic type. The fact that along the characteristic a special condition is valid 
shows that characteristics can pass information on the state of the fluid. The auxiliary variable, which has 
been introduced, is called celerity, and has a special meaning, it is the celerity with which a disturbance is 
propagated in water. Cunge et al. (1986) shows the relation between characteristic lines and the boundary 
and initial data requirements for the solution of Saint-Venant equations. In Figure 2.10 it is shown that 
in the case of subcritical flow at the upstream boundary of the computational domain one characteristic 
arrives from inside the domain, carrying information from another point of the computational domain. 
The other characteristic comes from the outside of the computational domain, hence it has to be replaced 
by a boundary condition. Such a boundary condition is given usually in the form of water depth, water 
velocity or discharge-stage relationship. The same analysis applies at the downstream boundary of 
the computational domain. In Figure 2.10b the case of supercritical flow is shown. It can be seen that 
the problem can be solved if two boundary conditions are given at the upstream boundary. The main 
conclusion is that for every characteristic entering a computational domain, one boundary condition has 
to be given. 
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(b) t (c) 
C1* 
LA cot R 
2 B.C. 1 B.C. 
LG. i LC. LC. ü 
Subcritical flow Supercritical flow Critical flow 


C* — pozitive characteristics 
Cc — negative characteristics 


Figure 2.10 Characteristic solution of Saint-Venant equations over a given computational domain; 
(a) Subcriticalflow conditions; (b) Supercriticalflow conditions; and (c) Critical flow conditions. 


Though this approach is very important, it cannot be applied to problems that include dissipative terms 
or for channels with irregular topographies. It is easier to use numerical approaches for finding the solution 
of Saint-Venant equations. 
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Chapter 3 
Discretization of the fluid flow domain 


3.1 DISCRETE SOLUTIONS OF EQUATIONS 


In the previous chapter the engineering problems were represented by equations. Many of the equations do 
not have an analytical solution. There are two approaches to find the values of the unknown function of the 
equations; one is to make further assumptions such that the equations are simplified and brought to a form 
which has an analytical solution; and another is to develop approximate solutions using numerical methods. 

The approach of simplifying the equations was successfully used in the past when computers were not 
available. The draw back of such an approach is that the solutions are valid just for a very limited number 
of simple engineering problems. 

The approach of using numerical methods has been used and developed extensively as computer power 
grew with time. In order to apply the method of finding a solution for a differential equation through numerical 
approximation, the computational domain (space and time) has to be transformed from a continuum domain 
into a set of subdomains, where the equations are still valid. The procedure of transforming the space and 
time of the computational domain into smaller parts is called discretization. With such an approach the 
infinite number of degrees of freedom of continuous equations is transformed into a finite number, which 
can be solved using a computer, because the continuous unknown variables are evaluated on the subdomain 
as discrete values. Hence the smaller the discretization, the smaller is the approximation error. 

There are different methods to discretize a computational domain, and in general the choice of the 
discretization method depends on the numerical method used for the solution of the equations. Each 
discretization method has its advantages and disadvantages. 

The discretization of the computational domain in smaller parts, is done by a so-called computational 
grid, defined by nodes (or points). Nodes are defining cells. The discretization of the space domain implies 
that the governing equations of the phenomena are written in a discrete form, that is, a linear system of 
algebraic equations is obtained that would give the values of the unknown variables at grid nodes. The 
linear system of equation is solved simultaneously (with few exception when this is not needed, as is the 
case of the explicit approach). 

In the process of transformation of the differential equations into a finite set of equations, through 
discretization, three stages are covered; the discretization of the definition domain of the differential 
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equation; the discretization of the unknown function in values at discrete points (i.e., the approximate 
representation of the solution); and the discretization of the equation itself. Each of these stages is described 
separately in this book. Discretization of the computational domain is covered in the current chapter, 
while the discretization of the solution and of the equation are explained in the chapters about finite 
differences and finite volume methods. Regardless of the numerical method used, different types of grids 
are used in engineering problem practice. Figure 3.1 shows an example of a simple grid, representing the 
discretization of Yongdam reservoir, in South Korea. The gridlines cross each other orthogonally, hence 
the grid is called Cartesian. 


Figure 3.1 Vertical Cartesian computational grid, the example of Yongdam reservoir, South Korea. 


The transformation of the continuous equations into discrete forms is done through methods such as 
finite difference, finite volume, finite element, boundary element, spectral methods, collocation methods, 
and so on. The first two are the methods used in this book and the presentation of grid generation is done 
in relation with these two methods. Discretization is done for both time and space, however the space 
discretization using in grids is more complex than that of time discretization, especially because it involves 
three independent variables (see Figure 3.2). 

Initially space grids were regular and structured, however as more complex numerical methods were 
developed, due to the computational power of computers, irregular physical geometries could be accounted 
for. Space grids, also known as mesh are of three types; structured, unstructured and meshless. Meshless 
discretization is not used by the finite differences and finite volume methods therefore not addressed here. 
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(a) Space domain discretization 


(b) Time domain discretization 


tinitial tfinal 


Figure 3.2 Conceptual discretization of the space and time domain. 
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Figure 3.3 Example of structured and unstructured mesh. 


The finite difference method requires a structured grid, while for the finite volume method both structured 
and unstructured meshes can be used (see Figure 3.3). 

For each cell of a computational grid a space is allocated in the computer's memory, because the 
unknown function is defined at the nodes of the grid or at the cell points of the grid. In case the grid 
extends over domains that are not of interest, the required computational memory is very high. 
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Several author presents comprehensive reviews on how meshes are generated and for which type of 
problems they should be used (Buell & Bush, 1973; Thacker, 1980; Ho-Le, 1988). 


3.2 SPACE DISCRETIZATION 

Space discretization is discussed for two and three dimensions of space, because the one-dimensional space 
does not poses any challenge in discretization. One-dimensional computational domains are discretized in 
a computational grid by selecting nodes unequally or equally spaced (see Figure 3.4). 


m | | =] >. 
h + + h 

Ax, AX» Ax Axn 

unequally spaced grid 
(b) 

| | | } | } } ». 
Jj. —— + I—4 

Ax Ax Ax Ax Ax 


equally spaced grid 


Figure 3.4 Conceptual discretization of the space and time domain. 


3.2.1 Structured grids 


Structured grids are simple, used for the discretization of a rectangular domain. Each of the axes of 
coordinates is divided by nodes in subdomains, independently of the number of subdivision on the other 
axes. For example in case of a 2D grid there are a number of nx and ny cells in the x and y direction, 
respectively. Each node or cell in the computational grid is identified by two indices. 

Structured grids are classified as: 


* Uniform regular grid, where the division of space in each direction is a constant (equally spaced) 
(see Figure 3.53); 

e Non-Uniform regular grid, where the division of space in each direction is not constant (see Figure 3.5b); 

e Orthogonal grids, where the grid lines are orthogonal to each other (see Figure 3.5a and 3.5b); 

e Curvilinear grids, where the grid lines are curves in space (see Figure 3.5c). 


Regular structured grids are especially used with the finite difference method, because the derivatives 
of the unknown function are easily computed with respect to all spatial directions. Moreover these types 
of grids do not require a transformation of the coordinates, as is the case with the finite volume method, 
or other numerical methods. The use of such grids is easily implemented in a computer programme, 
because the nodes and cells are determined using their indices, hence the structure of the solution matrix 
is known in advance. The main disadvantage of structured grids is that they cannot represent irregular 
complex domains. 
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Figure 3.5 Examples of structured grids. 


The only way for structured grids to partially cover irregular domains is by using the so-called 
multiblock grids, which are a combination of several structured grid blocks (see Figure 3.6). The solution 
computation is done on each block separately, however the computation advance from one block to another 
requires special treatment. 


Block 2 


Block 1 


Block 3 


Figure 3.6 Example of a multi-block grid. 


Apart from the multi-block type grids, which try to overcome irregular domains, there is another 
particular type of discretization grid called a staggered grid. These type of grids are specific to the problem 
that is solved by the differential equation and the way derivatives of the unknown function are expressed. 
They are used in case of a phenomenon that requires solutions for a system of differential equations and 
the different unknown variables have different approximation representations for the derivatives. For such 
grids different unknown variables are determined in different nodes of the grids. 

Anexample of a staggered grid in one dimension is shown in Figure 3.7 in which Saint Venant equations 
are solved for a one-dimensional flow in open channels. 


Downloaded from https://iwaponline.com/ebooks/book-pdf/650788/wio9781780400457.pdf 
bv IWA Publishina. publications@iwap.co.uk 


38 Computational Hydraulics 


tA e -unknown value 
o -known value 
' bo g & oq. 
n+l q 3 
n 
] ] ] i ] A 
L x 


Figure 3.7 Example of a staggered grid. 


3.2.2 Unstructured grids 


Unstructured grids are complex and are not generated easily, however they are highly flexible in 
representing the computational domain. The number of nodes and cells of an unstructured grid are not 
related by an equation, hence it is not easy to identify a cell using indices only, as in the case of structured 
grids. Each cell of an unstructured grid has to be defined by the nodes that are forming it, therefore all 
coordinates of the nodes have to be defined along with the information regarding the cell to which they 
belong. Unstructured grids are used for the finite volume method and they are the grid of preference for the 
finite element method. They have the advantage that they can represent some areas of the computational 
domain in a very detailed manner, however they have the disadvantage of longer computational time. An 
example of an unstructured grid is shown in Figure 3.8, as general representation. Figure 3.9 shows the 
structured and unstructured grid representation for a portion of Yellow River, in China. 


Figure 3.8 Example of unstructured grid. 
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Figure 3.9 Structured (a) and unstructured (b) grid representation for the middle section of the Yellow 
River, China. 


3.2.3 Grid generation 


Grid generation is a difficult task, and can be performed manually just for simple domains that require a 
small number of cells. For complex computational and/or large computational domains tools are used to 
create the grid of a given domain. Such tools are called grid generators. Grid generation is a field of study 
on its own and several methods are available to develop grids. These methods use algebraic relations 
to create the nodes and cells of a grid. Grid generators can take into account several requirements to a 
grid such as maximization of the cell area or orthogonality. In general grid generators are dependent on 
the type of differential equation to be solved, that is, the grid generators for hyperbolic equations are 
different from the elliptic ones, because of the type of boundary conditions that need to be defined at the 
nodes of the grid. 

In general grid generators produce a computational domain covered by a grid of triangles, by first 
creating all the nodes of the grid and then connecting them to form triangles. The best way to form these 
triangles is defined by the Delaunay triangulation, which ensures that very thin triangles are avoided 
(Boubez et al., 1986); or by using the Voronoi diagram (Brostow & Dussault, 1978). The Voronoi diagrams 
consist in building N polygons V(i) out of N given points, in such a way that each polygon is centered in 
point i, and the locus of the points most closest to point i is included in V(i). 
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The output of a grid generation is a set of elements along with a set of points that define it. The set of 
elements represents the grid topology. Depending on the approach taken to create the elements of a grid, 
the grid generators are (see Figure 3.10): 


* Creating first nodes then elements; 

* Creating first grid topology then nodes; 

* Creating nodes and elements simultaneously; and 
* Adapting grid from pre-defined grid templates. 


Topology Description 
Nodes creation gi 


Node description 


Grid topology — Grid smoothing 
creation 
Mesh generation 
methods m Grid based 
Nodes and grid Mapping 
topology 
Conformal mapping 
Use of predefined i Y 
templates — Geometry decomposition 


Figure 3.10 Classification of grid generation methods. 


From the four methods listed above the templates are generated separately, independently of the 
computational domain and they are later adapted for the given computational domain. The approaches 
known to create templates are: grid-based, mapped elements, conformal mapping and nodes and elements 
created simultaneously. 

There are several quality checks that need to be performed over a grid after it has been generated using 
a grid generator. The most important aspects are: grid smoothing, grid density and grid conformity. In case 
these three qualities of a grid are not fulfilled then the grid has to be refined manually. 

Grid smoothing refers to the shape of the elements. Smoothing is performed by repositioning the nodes. 
Repositioning is a difficult task that requires iterations, as repositioning might be wrong. 

Density of a grid refers to the ability of a grid to change from one part of the computational domain to 
another depending on the computational needs. Conformity is related to the density of a grid, because when 
a grid is refined smaller elements are introduced and they might be distorted. In general orthogonality is 
the measure to check for grid conformity. 


3.2.4 Physical aspects of space discretization 


There are problems for which the domain over which the solution is sought is bounded by special 
elements, such as structures, or may have a special form in the physical plane. Such special domains with 
special boundaries may be discretized to obtain the computational domain, which may be represented 
in a simpler format through a simple transformation of coordinates. After the solution is computed all 
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values of the unknown function are placed back to the initial computational domain, which represents 
the real domain. Such an example is a distorted L-shape region that can be transformed into a simple 


I-shape (see Figure 3.11). 
b n 
— 
X 


Figure 3.11 Computational domain transformations from one system of coordinates to another. 


Transformation should be done with care because it may introduce accuracy problems in the solution. 
In case of a bump a transformation grid may eliminate the bump (see Figure 3.12). 


/ N | 


Figure 3.12 Computational domain transformation for a bump. 


Fletcher (1991) emphasizes the importance of carefully making such transformations because the 
computational domain is distorted and so is the corresponding approximate solution. 


3.3 TIME DISCRETIZATION 


Discretization of time is separated from the one in space. Time discretization is used for problems, where 
the boundary of the computational domain varies in time, as is the case with free surface flows. In each 
time step values of the unknown function are computed at all nodes of the grid in space. When all values 
of the unknown function have been determined at time step n computation advances to the next time step 
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(n + 1). Numerical methods algorithms are written with respect to time level n and (n + 1), considering that 
all values of the unknown function are known at time step n, for all nodes of the space grid, and a new 
computation starts to determine all the values at time level (n + 1) where values of the unknown function 
are considered unknown (see Figure 3.13). 


t^ 


EOT 4 


0 1 2 3 kl i i+] 


e -unknown value at time level (n+ 1) 
o -known value at time level n 


Figure 3.13 Time discretization. 


Depending on the numerical method to determine the approximate solution of a differential equation, 
there are two main types of discretization with respect to time; explicit and implicit methods. 

In explicit discretization the values at the grid nodes of each unknown variable at the new time step 
of computation depend only on known values at the previous time step. The advantage of this kind of 
discretization is that in many cases no system of equation is formed, hence no solver is required, because 
the solution is straightforward computed step by step. Such a discretization is easy to implement in a 
computer programme. The disadvantage of such an approach is that solution becomes unstable if the time 
intervals for computation are selected bigger than a threshold value. Such a threshold value is usually 
determined from stability conditions of the solution combined with physical characteristics of the variable 
that needs to be determined. 

The implicit method of time discretization determines the values of the unknown function at the grid 
nodes from values of the function at neighbouring points at the same time level. An algebraic system of 
equations is formed, that requires a solver to be solved. The advantage of such method is that they are 
stable and no threshold value on the time interval of computation is required. From an accuracy point of 
view one should always try to solve the problems with small time intervals. 

In order to take advantage of both implicit and explicit methods of discretization for time, combined 
methods are also used. They are called semi-implicit methods, where the variable at the new time step 
depends both on the old and new time step. 
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Chapter 4 
Finite difference method 


4.1 GENERAL CONCEPTS 


The finite difference method (FDM) dominated computational science since its beginnings. This 
was the method of choice in computational hydraulics in the 1960s and 1970s. While other methods, 
such as the finite element method and the boundary element method have been extensively used, the 
numerical schemes developed using the finite difference approach are still used for a wide range of 
engineering problems. 

The finite difference method uses the differential form of equations representing a physical phenomenon, 
and solves the system of equations by determining values of the unknown variable function in discrete 
nodes of the definition domain. The domain is discretised in grids, as shown in Chapter 3. The values of 
the unknown variable are determined in the nodes of the grid. The value of the variable in any point of the 
computational domain, that is not a discretization node, is obtained by interpolation between the known 
values at discretization nodes. 

The finite difference method can be characterized as follows: 


* Utilizes spaced grids, mainly uniform ones; 

* At each node, each derivative is approximated by an algebraic expression that relates the values of 
the unknown function of the variable with reference to the adjacent nodes; 

* A system of algebraic equations is obtained and solved for the unknown variable, which is also 
known as dependent variable. 


The applied finite difference method is expected to generate estimates of the unknown function that 
will reproduce the analytical shape of the same function with sufficient accuracy. The accuracy depends 
on the selected nodes (points) of discretization and the calculation method of the unknown function 
derivative. 

In case time is one of the independent variable of the equation to be solved, three kinds of numerical 
schemes are used in practice: 


* Schemes in which the value of the unknown variable at time level n + 1 can be computed directly (or 
explicitly) from the value at time level n. Such schemes are called explicit schemes. 
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* Schemes in which the determination of the derivative involves the value of the unknown variable 
at time level n + 1. In such schemes there is an implicit relationship between the derivative and the 
unknown variable, hence they are called implicit schemes. 

* Schemes in which the explicit and implicit approaches are combined. 


Ordinary differential equations, depending on the computational approach of the numerical scheme, 
are classified as one-step methods and multi-step methods. If the unknown function at a node of the 
grid can be computed using information only from the previous node of the discretization grid then the 
method is called one-step. One-step methods are also known as self starting methods, because the given 
initial condition is sufficient to perform the computation for the value of the unknown function in all the 
nodes of the discretization grid. If the value of the unknown function in the current computational step 
is computed based on information from more than one preceding step of computation then the method 
is called multi-step. 

Each type of scheme has very different properties, as shown in Chapter 6 of the book. There are many 
finite difference methods available for solving differential equations, however due to the properties of the 
equations used in the aquatic environment, only few of them are applicable, therefore only a selection of 
these methods are presented in this chapter. 


4.2 APPROXIMATION OF THE FIRST ORDER DERRIVATIVE 


Given a function f(x) with one independent variable x, the definition of its differential with respect to x, as 
given in elementary calculus, is (see Figure 4.1): 


"too gm FO A9 - fo) 
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i i >X 
x x+Ax 


Figure 4.1 Definition of the differential. 


Nowadays calculations are carried out on computers and these cannot deal with a limit for which 
Ax — 0, therefore a discrete form of the continuous case is introduced. 

The function f(x) is defined by a set of values available at a number of discrete points of the domain of 
definition of the function. Figure 4.2 shows a set of discrete points x; where the function is known. The 
nodes x; are positioned at distance intervals Ax; = x; — x; The notation f; is used to denote the value of the 
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function f(x) at the i-th node, that is, f; = f(x). The set of nodes is called the computational grid (as already 
defined in Chapter 3). 


fox) 
fex) 


NE bod m 
XQ X2 m Xi Xie x 


Ax, AX AX 
Figure 4.2 Discretization nodes for a one-dimensional function. 


The approximation of the function's derivative is done using the discrete values of function fin a number 
of points of the domain of the function. This approximation is obtained using Taylor series expansions. 

Consider, for simplicity, an equally spaced grid of the independent variable x. A Taylor series is an 
expansion of f(x + Ax) or f(x — Ax) for a finite Ax at a fixed node x and can be expressed as: 


2 

f(x Ax) = f(x)  Axf'(x) + É f(x) + ar f")... (4.2) 
2 3 

f(x — Ax) = f(x) — Axf'() + a f(x) ie f" (x)... (4.3) 


where f" is the first derivative of f with respect to x, f”, the second derivative and so on. 
Rearranging (4.2) gives: 
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jav (x) + f(x) + «| = f'Q) + O(Ax) (4.4) 


If equations (4.2) and (4.3) are subtracted from each other it yields: 


fox Ax)- f(x- Ax). 1 
2Ax ~ 2Ax 


3 
parow * 2^7 pa) * «| = f'(x) + O(Ax?) (4.5) 


If Ax is small then it can be considered that: 


f(x) REF f(x + A) TU f(x) (4.6) 
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or 
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Equations (4.6)—(4.8) are examples of three simple finite difference approximations, as represented in 
Figure 4.3; the forward, backward and centered difference, respectively. 
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Backward derivative 
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Figure 4.3 Forward, backward and central approximations of the derivatives. 


In Figure 4.3 function f(x) is represented for a small interval near the discretization point x;. It can be 
seen that forward differencing uses the forward values of function f(x), that is, the values of the function 
in point (x, ,, f(x; + Ax,)), backward values in point (x; ;, f(x; — Ax,_,) for backward differencing and values 
on both sides of the function, in points (x;, f(x; + Ax,)) and (x, ,, f(x; — Ax; ;)) for central differencing. 

The term O(Ax) in equation (4.4) or O(Ax?) in equation (4.5) is a truncation of the Taylor series after 
the second term (right hand side of the equation). These terms are referred to as the truncation error (T.E.) 
because they represent the error when equations (4.6)—(4.8) are used. Hence the T.E. can be defined as the 
difference between the partial derivative and its finite difference representation. 

For smooth functions (the ones that have continuous higher order derivatives), and for a sufficiently 
small Ax, the first term of the truncation error, can give the order of magnitude of the error. Depending on 
the order of the term O(Ax), the following terminology is used: 


* Equation (4.3) is first order accurate in Ax 

* Equation (4.4) is second order accurate in Ax 

* Equation (4.6) is a higher order approximation to f'(x) than (4.3) or (4.4), which shows that if Ax > 0, 
(4.5) > f"(x) in a faster way than (4.4). It is important to note that the term is higher order not better 
approximation. 


Several higher order approximations of the first derivative of a function can be obtained by linearly 
combining the two Taylor series representations given by equations (4.2) and (4.3). 
In case that the discretization grid is not equally spaced (see Figure 4.2), equations (4.2) and (4.3) are: 


(Ax) 


3 Jd") dose (4.9) 


f(x; + Ax) = f(x) + Ax, f’ x) + —— a) + 
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(Ax 


fO = Ax) = fi) Aaf a) ++ SED” ra - 0 


A f” te (4.10) 


If equation (4.9) multiplied by Ax; ,, is added to equation (6.10) multiplied by Ax,, the following equation 
is obtained: 


1 
Ax. a AX; 4 fas. 
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Ax. —s 4G te 
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Ax; + Ax, Ax; f | Fx yo 


It can be noticed that in case of an equally spaced grid equation (4.11) is simplified to equation (4.6), 
that is, the cantered difference approximation is obtained. 

From equation (4.11), it can be noticed, that three values of the unknown function in three discretization 
nodes are used to compute the discrete value of the centered derivative at the node x,, that is, the nodes x; 
x; and x;,,. The truncation error O(Ax;Ax; 4) is of second order. 

Equation (4.11) can be used to determine higher order approximations for the first derivative of a 
function. For an equally spaced grid, the following expression holds: 


8(f(x; + Ax) - f(x; — Av) _ 


, e 
12Ax I (%) - 


a 


— f (x)- 5L f?(x.. (4.12) 


If equation (4.11) multiplied by 4 is subtracted from equation (4.11), written for a 2Ax grid spacing, 
equation (4.12) is obtained, which is a fourth order approximation of f(x). 

It is important to be noted that the validity of the Taylor series regarding the existence of higher order 
approximations of the derivatives depends on the existence of the derivatives themselves. If derivatives do 
not exist, then the higher order approximations do not hold. 

In a similar way the first derivative of a function of one variable is approximated by a finite difference, 
a first order partial derivative of a function depending on two or more variables can be approximated 
by finite differences. 


Figure 4.4 Representation of function g(x,t) in the (x,f) computational domain. 


Let g = g(x, f), than the first order partial derivatives of the function with respect to x or t are (see Figure 4.4): 
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Og _ &£u- 8l 
mo Me (4.13) 
and 
dg gi = gi 

NM. 4.14 
Ot "m At Pug 


where the following notations were made: 


dg 
Ox 
independent variables x and t, respectively; 

— gi is the value of the function g(x, f) at node (x; t"); 

— subscript indices are used for space, while superscript indices are used for time. 


"T and 7f an are the partial derivatives of function g(x, f) at node (x;f"), with respect to 


4.3 APPROXIMATION OF HIGHER ORDER DERRIVATIVES 


Differential equations of phenomena in water related sciences contain higher order derivatives. These are 
expressed using different approaches, such as polynomial fitting or Taylor series expansion, which provide 
a systematic way to estimate higher order derivatives of any order. It is important to note that the function 
for which the derivatives are approximated should be smooth. 

The approximation of higher order derivatives using Taylor series expansion is demonstrated here. The 
function is considered to vary with several independent variables, therefore partial derivatives are used to 
derive the expressions. The independent variable is considered to be x. 

Consider an equally spaced grid on which higher order derivatives are approximated in the node i (point 
of interest). A number of / points are considered at the left side of nodes i, and r at the right side of the node 
(see Figure 4.5). The value of function u in grid point i+ m (m=-l,...,7) based on value of the function in 
grid point i; and the higher order partial derivatives in the same point, is given by: 
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Figure 4.5 Equally spaced grid on the x axis. 


In equation (4.15), time level n has been fixed, hence the indices for the independent variables t are 
not shown, for simplicity of expression. 

If each of the r +/+ 1 equations of type (4.15) is multiplied by a constant coefficient a,, and then all 
are summed up: 
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$ r E r Aw r r (Ax)? 9?u 
> eel > «Je > man} at of 2 zt an) 2! ax]. 
m=—l, m#0 m--l, mz0 m--l, mz0 1 m--l, mz0 d 
r (Ax Qu Z (Ax) otu 
3 dion 2. 
| bj m «| 3! Qx? [S by mom 4! ox’ + 


m--l,mz0 m--l, mz0 i 


(4.16) 


Equation (4.16) is used to determine the r + / coefficients a,, of the derivative depending on the desired 
order. For example if the numerical solution requires a first order derivative, fourth order accurate the 
values of the coefficients b, are set as follows: 


b, -1 
b, =0 (4.17) 
b, = 0 
b, =0 
where 
b, = 3 m*a, (4.18) 
m--l,m 


The system (4.17) gives four equations that would require the use of four discretization points, in order 
to have a unique solution of the system. 

In general if the k-th derivative is needed then the highest derivative to be neglected is of order k + p — 1, 
that is, k + p — 1 grid points are needed to determine the solution, and equation (4.18) can be re-written as: 


b, = Oy = M oma qol. kep-l (4.19) 
m--l,m 
with 
ô, = l; fg =k (4.20) 
i 0 ifq#k 


with ô, known as the Kronecker delta. The solution exist and it is unique if: 
l+r=k+p (4.21) 


Based on (4.19)—(4.21) the order of the T.E. is determined by multiplying the next higher derivative in 
the truncation error series with the coefficient: 


Dyas = Y mPa, (4.22) 


m=-l,m 
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Based on (4.18), selecting the centred difference approximation, and considering an equally spaced 
grid, it can be noticed that the second order derivative of a function u(x, f), for which the domain is 
represented in Figure 4.6, can be computed as: 


2 n n n 
ou ul, — 2u” + 2u}, 


Ix hey == (Ax? (4.23) 


AX, AX 


At 


Figure 4.6 Domain of definition and grid discretization for function u(x, t). 


Another way of developing finite difference approximations are, as previously mentioned, polynomial 
fiting, such as splines or Lagrange interpolations. These approaches consist in fitting a polynomial of 
a certain degree to a set of selected nodes of the grid. For example, in case of a function u(x), of one 
independent variable, the slope of the function u, at node x; is approximated by considering the derivative 
of the considered polynomial at the point (x;, uj). The development of these approaches are presented in 
detail in Boyd (1989), Butcher (1987) and Duran (1999). 


4.4 FINITE DIFFERENCES FOR ORDINARY DIFFERENTIAL EQUATIONS 
4.4.1 Problem position 


Ordinary differential equations (ODEs) contain a function that varies with respect to one independent 
variable only. They contain derivatives of the unknown function and/or the function itself. The general 
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form of an ODE is: 
H2 du E (4.24) 
w = n o fut) 


where u = u(f) is the unknown function, t is the independent variable and f(u, f) is an arbitrary function of 
u(t) and t. 
Integrating equation (4.24) yields: 


u(t) = [re t)dt 4 C (4.25) 


where C is a constant of integration. There are an infinite number of values for C, hence the obtained 
solution of equation (4.24) is not unique. However in the daily practice of engineering solving problem, 
there is interest just for a unique solution for a specified range of t, over a given interval. In order to choose 
a solution an extra condition, called initial condition, must be specified. 

The problem is formulated as follows: Determine the unknown function u(f), over the interval (r'"'i4!, 
tn^) (see Figure 4.7), which satisfies the following set of equations (the given ODE and the initial condition): 


Pi (4.26) 


u( tia! ) = yinitial 


u(t) 
U(E) |j- M N 
4 LJ 
ut) Pee d 
ra ' - 
f utes 
, nei i 
2^ u i 
L H 
«Ue T 
ue j 
i f >t 
(initial. t^ tn? (final 
— 
At 


N intervals 


Figure 4.7 Discretization domain and solution approximation for function u(t). 


The problem formulated by (4.26) is called an initial value problem (IVP), for the ODE, and has a 
unique solution (Asher & Petzold, 1998). 

Examples of engineering problems in water related field that are represented by equations like the ones 
of (4.26) are the steady flow in open channels; unit hydrograph or reservoir storage equation. 

There are many more numerical schemes available for ODEs, such as explicit, implicit, Adams— 
Bashford, or Adams-Moulton (Ascher & Petzold, 1998). A selection of the most used numerical solutions 
for ODE, as they are used in water related problems, is presented in this section of the book. 
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4.4.2 Explicit schemes (Euler method) 


The problem defined by (4.26) can be solved using a discretization of the independent variable (i.e., time) 
between the values ft’ and f/7«/; The discretization is made into N intervals of width Aż (see Figure 4.7). 


ro [initial +n- At (4.27) 


where At = (tne! — initia JN: n=0, 1,...,N; and u(t") =u" is the value of u(f) at the n-th time-step (see 
Figure 4.7). 
In order to approximate the derivative du/dt, the forward difference expression for the derivative 


™ t TÅ aa DEA) (4.28) 


is substituted into equation (4.26). It yields: 


u”+! — u” 
~ — + O(At) = flu", t") (4.29) 
At 
Re-arranging expression (4.29) gives: 
yr! zu + At. f(u”, t”) (4.30) 
It is important to note that O(Ar) term has been neglected and by that the equation is now an 
approximation. Equation (4.30) provides an expression to calculate u(f) at a particular time-step r’*! from 
its value at a previous time-step, t". This is known as explicit differencing numerical scheme. 
If u? = uii! is known than u! can be estimated; from u!, estimation of u? can be done, and so on. After 
N steps, u” = uf"! is determined. 
In each approximation of the derivative term O(A?) is neglected therefore there is a small error, that is 
accumulating in time, hence, the error in the calculation of u! is carried over in the calculation of u*, and 
so on (Figure 4.8). 


Kit) + 
E a “fun” 
jem fin 
i f >t 
t? t! t? ee t^ 


Figure 4.8 Error accumulation. 
The step by step procedure of computing the values of u” in time can be expressed as: 
uy" ad u? 3 AL x Siw, t') (4.31) 
j=0 


The explicit method approach for solving the ODE (4.26) is a one-step method approach. 
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4.4.3 Implicit schemes (Improved Euler method) 
Using the same notations as in the previous section, in order to derive the implicit scheme, the term du/dt 
term in equation (4.26) is replaced by a forward differenced approximation and the value of the function 
flu, À is computed at the n + 1 node of the discretization. This yields: 


n+l 


"MA C O(Ar) = f, in) (22) 


Re-arranging expression (4.32) gives: 
ut "n At . f(u’, t") = u” (4.33) 


This formula allows to calculate the unknown value of u” at a certain time-step from its value at next 
time-step, u”*!. Unfortunately to do this, value of u”™! needs to be known and replaced into the function fin 
the first step of computation. This is what is known as an implicit differencing scheme and usually cannot 
be solved within a single step when the expression of the function fis complicated. Seldom it is necessary 
to use a technique called iteration, such as Newton's method. 


4.4.4 Mixed schemes 


If the explicit and explicit schemes are averaged together, the first-order ODE, defined by (4.26) is 
differenced as: 


um — y” + O(AD = lrae i") + fu! p] (4.34) 
At 2 i > 


By rearranging expression (4.34), it yields: 


unt -- 


m u" -— sli t”) + fu”, [pH ] (4.35) 


Expression (4.35) is more complicated than the explicit or implicit schemes and may require iteration 
to be solved. It is known as trapezoidal differencing scheme or a Crank-Nicholson scheme. 


4.4.5 Weighted averaged schemes 


A more general formulation of the explicit, implicit and mixed schemes is the weighted averaged formulation, 
which gives one single formula for all three formulations. Consider as point of interest £, which is not 
positioned at a discretization point. Position of 7? inside the interval (r", t") is defined by the weighting 
parameter 0, which gives its position with respect to the nodes of the discretization grid (see Figure 4.9). 


t? — ff 
- 4.36 
? At (4.36) 
with 0 ranging from 0 to 1. 

Approximation of equation (4.26), using point P gives: 
un! — u” up 
LO V m t 4.37 

apo 8 Sul?) (4.37) 
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tni " " nei >t 
9At | (I-9Atj 


At | 


Figure 4.9 Weighted average position of a Point P in the interval (t, t”). 


The value of f(u, À in point P is calculated by linear interpolation between the end points of the interval as: 


fu, t?) = 0— 0)f(u", 1) + 0- fur) (4.38) 


which gives the final expression for u”*! as: 


ut! = u” + At - ü = 0) i flu", t”) +At-0- fu, poH) (4.39) 
u(t) 
u( tnt) 
ut! 
a ; li At “fe ) 
u 
> 
t^ t^? 
At 
fut fwd fa^ 
fart ppt 5 fant p 
fasi) | fait”) 
> | » » 
t^ tn t^ t^^ t^ toH 
At At At 


Figure 4.10 Integration of function f(u, t) over interval At. 
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It can be noticed that scheme defined by expression (4.39) is: 


e explicit Euler scheme, if 0 = 0; 
* implicit scheme, if 0— 1; 
* mixed scheme, if 0 — 0.5; 


Coefficient 0 can have any value in the interval [0, 1], hence many more schemes can be defined. One 
of such scheme, known as Galerkin scheme is obtained for 0 = 2/3. 

If the variation of function f(u, t, from equation (4.26), is represented graphically, for a given u(f), then 
it can be noticed that the difference between the three numerical schemes formulations, explicit, implicit 
and trapezoidal, is in the way the integration of f(u, f) is done over the interval At. From Figure 4.10 it can be 
concluded that explicit approach is underestimating and implicit approach is overestimating the value of f(u, f). 


4.4.6 Runge-Kutta methods 


A similar approach to the weighted averaged schemes are the so called Runge-Kutta methods in which a 
set of multiple internal points of computation are selected inside the interval (¢”, t+). 

The general numerical scheme for Runge-Kutta methods, for a number of m points inside the 
computational domain, can be formulated as: 


k, = fus") (4.40) 


j-1 
k = e tA: Saiki I" 3e 


l=1 


where m 2 1, and w, c;, and a; are numerical coefficients. Depending of the number of internal points m, 
the corresponding Runge—Kutta method is called m-stage. 

Legras (1971) presents formulas for Runge-Kutta coefficients, up to 8-stage. One of the most popular 
Runge-Kutta formula is the one in 4-stages for which the formulas defined by equations (4.40) are: 


kı = f(u", t") 
k, = (e thy xj 
"m " (4.41) 
t 
k, = (e thon 3 
k, = f(u" + k At, t” + Ar) 
with the solution at time step n + 1: 
+1 At 
ue =u" + * (k, + 2k, + 2k, + k4) (4.42) 


There are many more numerical schemes available for ODEs, such as Adams—Bashford or Adams— 
Moulton (Ascher & Petzold, 1998). 
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4.5 NUMERICAL SCHEMES FOR PARTIAL DIFFERENTIAL EQUATIONS 


The principle of the finite difference methods is the same for PDEs as for ODEs, however PDEs are 
more complex to solve than the ODEs, because the unknown function that needs to be determined is 
differentiated with respect to both time and space. The boundaries, and the conditions which the solutions 
must satisfy on the boundaries, have considerable influence on the solution of a PDE. The influence of 
boundary conditions on the solution of PDEs is very important because a PDE might have solutions for a 
certain set of boundary conditions and might be insolvable for another set. 

As it has been shown in Chapter 3, PDE's are of three kind. Each of them is be treated separatly 
in this section of the book. Due to the special importance of the hyperbolic and parabolic equations 
in fluid dynamics, emphasis is given to these two kind of equations. Principle of the FDM approach is 
demonstrated on 1D types of equations, which is the same for 2D and 3D. In case of PDEs, 1D refers 
strictly to the space component, time independent variable is the second independent variable with respect 
to which the unknown function is defined. 


4.5.1 Principle of FDM for PDEs 


The solution of a PDE cannot be determined analytically, if the PDE is defined on a finite domain. Numerical 
methods determine a solution of a PDE by replacing the unknown function continuous derivatives with 
discrete approximations. The solution is found at a number of discretization grid points in space and time 
(see Figure 4.6). The numerical solution of a function u(x,f), varying in one dimensional space and in time, 
is determined at grid points x, = iAx, and time points nAt. The notation mentioned in equations (4.13)— 
(4.14) is used throughout the chapter, that is, u? = u(iAx, nAt) = u(x;, t"). 

Any approximation that is used, in order to determine a solution for the engineering problems should 
be a good one hence it must be consistent, stable, and convergent to the analytical solution. This issue is 
detailed in Chapter 6 of this book. 

The principle of the FDM for PDEs is demonstrated considering the simple PDE of the scalar linear 
advection equation, which is of hyperbolic type: 


P uc ed (4.43) 


where a is the advection speed, assumed constant and positive. Appropriate boundary conditions are assumed 
to be specified at the upstream boundary. The unknown function u(x, f) is one dimensional in space. The 
equation is linear because it consists of a sum of separate terms involving the advected quantity, u. 

When the advection speed a is a constant, it can be demonstrated that any initial value of u is advected 
unchanged along the x axis. The general solution of the equation (4.43) is written as: 


u(x, t) = F(x — at) (4.44) 


where F(s) is an arbitrary single-valued differentiable function; with s = x — at. Substitution of expression 
(4.44) into the left-hand-side of equation (4.43) yields: 


— lg "P3. 


OF(s) " OF(s) _ dF ds " dF ðs | dF( Os Os 
Ot 4 Ux dst "^dsóx ds 


]-z ata+l=0 (4.45) 
ds 


Deduction (4.45) shows that u(x, f) is carried along x at a constant velocity a, and the shape of u(x, f) 
remains unchanged (see Figure 4.11). 
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Figure 4.11 Advection of u(x, 0) along the x axis. 


In order to solve equation (4.43) using the discrete approach, the time-derivative between points (i, n) 
and (i, n + 1) are used (see Figure 4.12): 


= — 4.46 

Ot At" pes 
along with the space-derivative: 

Qu _ uj =u (447) 


Ox Ax, 


i 


Computational Domain 


AX; AX; 


B.C. - Boundary Conditions 
I.C. - Initial Conditions 


Figure 4.12 Unequally spaced grid discretization for the advection equation computational domain. 
Replacing (4.46) and (4.47) into (4.43) yields: 


ur = Cr” ut, + (1 — Cr") - u” (4.48) 
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(4.49) 


Expression (4.49) is called the Courant number, and depends on the grid size at time level (n + 1) and space 
size at (i + 1). Assuming that the values for u are known at any point i of the space, at time level n, the value 
for u at any point in space, at time level n + 1 can be computed using (4.48). Equation (4.48) cannot be applied 
at the left-hand boundary of the domain where i= 1, because there is no point to the left of the left boundary. 
This analysis shows the need to know the Boundary Condition, ut! , at the left boundary (see Figure 4.13). 


Computational Domain 


Known Left B.C. 


Given I.C. 
B.C. - Boundary Conditions 


I.C. - Initial Conditions 


Figure 4.13 Definition of boundary conditions for advection-equation. 


The notion of stencil is introduced here by definition. A stencil is the representation in the computational 
domain coordinates, of the nodes that relate to the point of computation by using a numerical approximation 
method. In case that two independent variables, space and time, are used to express a function in a 
differential equation, than in the notation used for stencils, n — 1, n, n + 1 indicate the time steps, and i — 1, 
i, i+ 1 indicate space step. It is always assumed that all values of the function under consideration are 
known at time steps n and n — 1, and the time step n + 1 is to be calculated. 

The stencil of the numerical scheme (4.48) is shown in Figure 4.14. 


or 


(i-1) (i) 


Figure 4.14 Example of a stencil for a numerical scheme. 
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4.5.2 Hyperbolic PDEs 


Demonstration of the application of the FDM for hyperbolic equations is shown on the linear advection 
equation. As seen from the principle of the method the first step in applying the finite difference technique 
is to discretise the problem domain. Next step is to find a suitable approximation to the hyperbolic equation 
for the considered discretization. 

The linear advection equation has two first order partial derivatives, du/dx and du/dt. There are 
many different numerical schemes to determine the appropriate solution of the equation, depending 
on the selected approximation of the two derivatives, with respect to time and space. If only the 
three different approximations for such derivatives (forward, backward and central differences) are 
considered a minimum of nine different approximations to equation (4.43) can be expressed. The 
main question is which approximation should be used, because while seeking a solution one needs to 
be sure that the differencing scheme should be as accurate as possible, numerically stable and easy to 
implement. 


4.5.2.1 Explicit and implicit schemes 


As already mentioned in the beginning of the chapter, a scheme where the value at the unknown time level 
(n + 1) can be computed directly from the previous time level is called an explicit scheme. Equation (4.48) 
demonstrated for the principle of the method for PDMs is an explicit type of scheme. Implicit schemes are 
further detailed bellow. 

Given the domain of definition of the unknown variable u(x, f), as defined in Figure 4.12, the first order 
space derivative of the unknown function is estimated using the values at time level n + 1, as: 


ntl ntl 
du uf — ut 


ox = AE (4.50) 
Substituting expression (4.50) and (4.46) into (4.43) yields: 
-Crhlum + (1+ Cg). utt! = u? (4.51) 


The stencil of the numerical scheme (4.51) is given in Figure 4.15. 


Figure 4.15 Stencil of the explicit scheme for advection equation. 


The unknown solution at point i for time level n + 1 depends on the known values at the previous time 
level n, as well as on the unknown value at point i — 1 at time level n + 1. The values of u at time level 
n+ l are linked together through a linear relationship. If equation (4.51) is written for all the nodes at 
time level n + 1, a system of equations is obtained and must be solved. The upstream boundary condition, 
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uj*! is known and, už, is known from the previous time step, hence it is possible to compute, u*!. Sweeping 


along i in ascending order provides, u7*', u7*', and so on (see Figure 4.16). 


* Known values at time level n 
x Unknown values at time level (n1) 
O Given initial and boundary conditions 


Stencil of 
the 
explicit 
scheme 


Figure 4.16 Explicit solution approach for hyperbolic equations. 


A scheme where there is a relationship between several points at the unknown time level, leads to a 
system of equations. The equations of this system must be solved simultaneously to compute the unknowns. 
Such a scheme is called an implicit scheme. If the relationship between the unknown points is complex (for 
instance, non-linear), the system must be solved in an iterative way. 

The differential approximations of the derivatives that are forward, backward and centred space are 
abbreviated as FS, BS and CS respectively. Similarly the ones for time are abbreviated as FT, BT and CT. 
Any numerical scheme that uses a combination of any of these two approximations is also referred to as 
a combination of the two abbreviations. The numerical scheme (4.49) is referred as FTBS, and (4.51) as 
BTBS. 


4.5.2.2 Forward time — centered space scheme (FTCS) 


The forward time, centred space scheme for the advection equation, based on the FT discretization of the 
time derivative of u, and CS discretization of the space derivative is: 


um — u’ ul up, 
miu o a-h 4.52 
N +a JA 0 ( ) 


or, re-arranging 


l aAt 
up = ur = 5 . A uh a5 ui) (4.53) 


Durran (1999) shows through Fourier analysis that this scheme is unconditionally unstable. Numerical 
models still use a FTCS scheme to calculate the few first time-step, because the errors introduced by using 
the FTCS scheme at the beginning of the computational steps, are negligible when compared with all the 
other errors that are present in a model. 
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4.5.2.3 Centered and upwind schemes 


The schemes define using the backward space approximation for the space derivative are called upwind 
schemes because the space information that is used to make the computations is coming from ‘upstream’ 
of the computational domain. 

The schemes that are using centered derivatives, either for space or time are called centered schemes. 
One of the most interesting schemes of this kind is the CTCS scheme: 

= up! Ui — Ui 


i . - 4.54 
2AI +a Ax 0 ( ) 


or, re-arranging 


At 
- PE Qd, — uta) (4.55) 


m =i ut 


u 


Scheme defined by relation (4.55) can produce two solutions to the same problem under special 
circumstances of boundary conditions. The stencil of the CTCS scheme is given in Figure 4.17. It can 
be seen that for a particular grid-point i and time-step n, the CTCS approximation equation (4.55) links 
together the values of U at (i, n + 1), (i — 1, n), (i + 1, n) and (i, n — 1). The value of u at grid point (i, n) is 
not part of the relationship. The same equation for the same grid-point i at the next time-step, n + 1, links 
(i, k+ 2), (i € 1, n * 1) and (i, n). Figure 4.17 shows the nodes for which the values of u effect the value of u 
at node (i, n) and those that effects the value of u at (i, n+ 1). 


t 


Figure 4.17 Relation between nodes for the CTCS numerical scheme. 


4.5.2.4 The Preissmann scheme 


In case of open channel flow modelling a frequently used numerical scheme is the Preissmann scheme, 
or ‘box’ scheme, or 4-point scheme. The scheme is implicit, very flexible and suitable for demonstrating 
some of the properties of the numerical solutions of hyperbolic types of equations (Abbott & Basco, 1989; 
Preismann, 1961). The scheme is developed based on four nodes of a discretization grid, as shown in (Figure 
4.18). The derivatives are not approximated at the nodes of the grid, but at points located inside the grid. 

Consider a grid equally spaced in space and time and the advection equation defined by (4.43). For a 
point P located between the points i and (i + 1) in space and time levels n and (n + 1) the following formulas 
are used to express the derivatives of an arbitrary function f(x, f): 


n 


ntl n n+l 
ou ur — uti unt! — ul 
t 


on, V Ar. CU CMT 


(4.56a) 
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B (4.56b) 


where y and 0 are weighting coefficients ranging from 0 to 1. 

Formulas (4.56a, b) represents a weighted average of the derivatives between points (i — 1, n), (i, n), 
(i, n+ 1) and (i — 1, n + 1). If equation (4.43) is valid for a flood wave travelling in an open channel, then 
(1 — y) is the distance between point P and the cross-section positioned at grid point (i — 1) (Cunge et al., 
1986). For 0 = 0, the scheme is explicit, and it becomes fully implicit for 0 = 1. This scheme is unstable 
for 0 « 0.5. 

Substitution of relations (4.56a, b) into (4.43), yields: 


+1 n 


n qn n+l _ „n+l 
HE de afa 8). 4 a pg a - 0 (4.57) 


ntl n 
Yin 7 Uia 


At 


ur — 
T( 0— : 
a-y) "Y 
If all values of u(x, f) at time level n are known then equation (4.57) relates the values of the unknown 
function u(x,t) at grid points (i, n + 1) and (i — 1, n + 1) as follows: 
ult! = a; -uti + b,- ut +c un (4.58) 


with coefficients a;, b;, c;, given by: 


' 1-wv-0-.Cr 
» 1-V-ü-8) Cr 
p= 0.C 
ld ibd (4.59) 

p- te Cr 
! ]-w-80:Cr 

At 
r= ay 


In case that a canal of length L is considered, over a time period T, and an equally spaced grid 
discretization of the computational domain in (M — 1) space intervals Ax, and (N — 1) time intervals Ar. 


Downloaded from https://iwaponline.com/ebooks/book-pdf/650788/wio9781780400457.pdf 
bv IWA Publishina. publications@iwap.co.uk 


Finite difference method 63 


Such a discretization results in M points along the x axis and N time levels of computation, as represented 
in Figure 4.19. Equation (4.57) is valid on every point i = 2,..., M, at a specific time level n. At point i = 1, 


the value of function u(x, f) is the upstream boundary condition and it is known, as are known all the values 
of u(x, 0) as initial condition to the problem. 


t 


Given B.C. 


Given B.C. 


Qh on 


| Given I.C. L 


Figure 4.19 Computational domain for Preissmann scheme. 


Writing equation (4.58) on every point i = 2,..., M, a system of (M — 1) equations, at time /"*! is obtained. 


Adding the known boundary condition u7*!, a linear system of M equations with M unknowns is required 
to be solved in order to have the solution. 


4.5.2.5 The Abbott-lonescu scheme for free-surface flow 


The second scheme that is often used for free surface flows is the Abbott-Ionescu scheme, which was 
developed to solve the system of one-dimensional non-conservative shallow water equations (Abbott, 1979): 


Ou Ooh g Oh 
hA ust (gh Ey 8 


(4.60a) 


oh Qu "T 
Ex cH Spe c E) M 


(4.60D) 
where g is gravitational acceleration; A channel water depth; and u flow velocity. 
The discretization grid used for solving these equations is a staggered one that is, the space- and time- 
derivatives of the variables A and u are not computed at the same grid point (see Figure 4.20). 


f 


Given B.C. 
Given B.C. 


Given I.C. L 


Figure 4.20 Computational domain for Abbott-lonescu scheme. 
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The derivatives of the equation (4.60a) are defined as: 


Oh hU — hp 
ot At 


du |] WD — uja Ui! Ufa 0 
At At 


Jk 2 2Ax T 2AX 


n+l _ „n+l n _ an 
du [5 Uii , Uia ao 


The derivatives of the equation (4.605) are defined as: 


n+l n 
Ou — wh ui 


ot At 


dh 1 hy Ey hi hi — hii, 
ài al cs ae 


dh _ 1 (kg-a | the He 
ox 2 2Ax 2Ax 


Figure 4.21 illustrates the principle of the Abbott-Ionescu scheme. 


(n+1)-9 


(11) OQ thee Qn E 
(i-1) i (41) (+2) 
PE ad ell 
i * i 
l EA i 
! Pd es 1 
PM i 
TO O i 
For continuity For momentum 
equation equation 


Figure 4.21 The Abbott-lonescu numerical scheme. 
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(4.61d) 


(4.6le) 


(4.61f) 
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4.5.3 Parabolic PDEs 


The solution of parabolic equations using FDM approach is demonstrated on the one-dimensional diffusion 
equation, which has the general form: 


Ou Q?u 
Lae lila 4.62 
Ot n oe) 


where D is referred as the diffusion coefficient. 

The same principle, as in case of hyperbolic equations apply, that is, a solution for the unknown function 
u(x, t) is determined at time level (n + 1), based on the values previously determined at time levels 1, 2,..., n. 

In the practice of solving parabolic equations two types of boundary conditions are common to be 
defined for the domain of computation; Dirichlet boundary conditions; and Neumann boundary conditions 
(Fletcher, 1991). The Dirichlet boundary conditions are when the unknown function is defined at both 
ends of the computational domain, while Neumann boundary conditions are when the spatial derivative of 
the unknown variable is known at the both ends of the computational domain. An initial condition is also 
required in order to find the solution of the diffusion equation (Figure 4.22). 


Computational Domain 


g g 
m m 
ae) gz 
D 9 
a .H 
=] z 
o o 
9 9 
4 4 
x 


Required I.C. 


B.C. - Boundary Conditions 
I.C. - Initial Conditions 


Figure 4.22 Computational domain, boundary and initial conditions for the diffusion equation. 


4.5.3.1 Explicit schemes 


The parabolic equation (4.62) involves a second-order space derivative. According to relation (4.23) at 
computational time level (n + 1), the values of the unknown function in three points in space are needed to 
estimate the second-order space derivative. This is a CS approximation. 

Consider an equally spaced computational domain (see Figure 4.23), a forward time approximation 
for the first order time derivative of u(x, f), and a central space approximation for the space derivative 
(equation (4.23)). Replacement of the approximate derivatives in (4.62) yields: 


n+l n n n n 
U; Uj ua — 2u; + Qu 


m Te. - (4.63) 


Re-arranging gives: 


ut =r- uty +(l—2r)- u; +r- uta (4.64) 
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(4.65) 


X 


[Ax Ax, 


Figure 4.23 Stencil of the explicit FTCS scheme for a parabolic equation. 


The stencil of the numerical scheme given by equation (4.64) is shown in Figure 4.24. 


Figure 4.24 Stencil of the implicit scheme for a parabolic equation. 


Other explicit schemes for the diffusion equation are the DuFort-Frankel or Three level scheme for 
which the corresponding equations are (4.66) and (4.67), respectively: 


2r 1-2r 
ntl _ n n n-l 
u -[3 Jes te [TE Je (4.66) 
ut = | EY us -| Y lu H "| (= Bur Bul) (4.67) 
i l+y J?^ (1+7)! ltr}? i i 


B--05-y-i (4.68) 


Ly uf = uja — 2u? tuj, (4.69) 
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The value of the coefficient yis arbitrary and has to be chosen with care so that in combination with the 
value of r fives a stable numerical scheme. 

Equations (4.66) and (4.67) involves values of the unknown function at three time levels (n — 1), n, and 
(n + 1), except for the case r = 0.5 when equation (4.66) reduces to the FTCS scheme. 


4.5.3.2 Implicit scheme 


The implicit scheme for diffusion equation is obtained applying the FT derivative at time level n, in point i 
of the space, and the CS discretization of the derivative applied at time level (n + 1) in point i of the space. 
The obtained FTCS implicit numerical scheme is: 


ru! -(ü-c-2rn-:u"* -r.uz-u (4.70) 


The stencil of the relation (4.70) is shown in Figure 4.24. 

The substitution of the derivatives in equation (4.62) leads to a system of equations, because numerical 
scheme (4.70) has to be written in all i nodes of the discretization grid, at time level (n + 1). The resulting 
system of equation, that provides the final solution is: 


d+2r) -=r e -| [ug] pe 
=r (1+ 2r) =r e| rur u^ 
=|... (4.71) 
-r A+2r) -r| uw | ns | 


Similar to the case of explicit schemes, a more general formulation of the implicit scheme is the three 
level scheme: 


(Ly): Qd = ut) — 7 (ut = uty) 
At 


a-[(— DL, -ur + BL,, -u]-20 (4.72) 


A particular Three level scheme, the one for which y=0.5 and B= 1, is often used for solution of 
problems. 


4.5.4 Elliptic PDEs 


As for previous types of PDEs, the application of the FDM method for elliptic type of PDEs is shown on 
the example of Laplace equation. Poisson equation is an example of an elliptic equation, and it is written 
in the form: 


Qu o 
35 * E = f x, y) (4.73) 


In case that the function f(u, x, y) = 0 the equation is called Laplace. 
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Notice that time is not involved in the elliptic type of equations, because they are for steady state 
problems, hence time independent. For the numerical solution the domain represented in Figure 4.25 is 
used and the partial derivatives are substituted by the finite difference approximation (4.23). Because 


time is not involved in this kind of problem, only subscript notations are used, the ones for space 
(i and j). 


(b) 


Required B.C. 


Required B.C. 
Required B.C. 


Required B.C. 


Figure 4.25 Domain of computation and stencil for elliptic type of equations. 


For an equally spaced grid, (i.e., Ax = Ay) ,the numerical scheme obtained after substitution is: 
Wn; Wu; + Hj T Hu; T Au; = 0 (4.74) 


The stencil corresponding to numerical scheme (4.74) is represented in Figure 4.25b. 

Similar to the implicit methods for parabolic equations relationship (4.74) needs to be written for all 
interior points of the domain and the result is a system of equations that gives the solution for the unknown 
function. Relation (4.74) is often referred to as the Laplacian difference equation. 

Analysing relation (4.72), it can be noticed that boundary condition along all the four edges of the 
computational domain must be specified in order to obtain a unique solution of the system of equation. 

Boundary conditions at which derivatives of the unknown function are specified makes the problem 
solving more complex. 


4.6 EXAMPLES 
4.6.1 ODE: Solution of the linear reservoir problem 


One of the most used ODE in hydraulic engineering is the linear reservoir. The linear reservoir is the continuity 
equation written for an inflow and outflow that changes the volume of water contained in a reservoir. The 
equation applies in reservoir operation, hydrological routing of flow, groundwater flow, and so on. 

Given a reservoir, as shown in Figure 4.26, the change in volume is given by: 


dV) _ 
dt — 


Q;,() - Q,, (t) (4.75) 


in which V(f) is the water volume stored in the reservoir, Q,,,(f) is the inflow into the reservoir and Q,,,,,(0) is 
the outflow from the reservoir. 
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Figure 4.26 Scheme of a linear reservoir. 


Equations (4.75) has two unknowns, the volume of the reservoir, V(f), and the outflow from the 
reservoir Q,,,,(f). Taking into account that the outflow from the reservoir is a function of the storage of 
the reservoir: 


Qon O = FV) (4.76) 


If the assumption is made that f(V(r)) is a power function then it can be written: 


Qu) = k- [VO] (4.77) 


where k is a storage constant, which shows how fast the storage is emptied and is measured in units of 
1/time; and a is a constant parameter. For a particular reservoir both a and b need to be determined from 
specific particularities of the reservoir. 

In general it is considered that a= 1 is a good value (Henderson, 1966; Sturm, 2001), in which case 
equation (4.77) can be rearranged as: 


VO) = 10,0 (478) 


If the simplified notation Q,,,(f) =Q and Q 
mass conservation equation (4.75) gives: 


(f) 2 I is used; and replacing equation (4.78) into the initial 


in 


ASA a7 (4.79) 
a 1? 

Re-arranging gives: 

1 dQ 

mu. NS M 4.80 

ka E (4.80) 

Or 

dQ 

UE = 4.81 

"SLE (4.81) 
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Equation (4.81) can be integrating applying the integrating factor method, for solving ODE. The 
equation (4.81) is multiplied with the term e^, hence: 


et a + ke*Q = ke*I (4.82) 


The right-hand-side of the equation (4.82) represents the derivative of the function Qe", hence integrating 
the equation, yields: 


eQ = frets dt (4.83) 
The final expression for the outflow is: 


D k |e*1 dt (4.84) 
e 
Expression (4.84) 1s not easy to be used in practice, because /(f) is a function that has discrete values. 
When /(f) is known the solution for Q(r) is determined by using numerical integration of the term e"7. 
In the particular case that (f) = 0, equation (4.81) reduces to a simple ODE, that can be numerically 
solved using finite difference methods. Assuming that Q(t = 0) = Q, the equation that is solved numerically 
is (see Figure 4.27): 


dg 
& = -kQ (4.85) 


T T T 
— Analytical solution 
45 | -- Implicit solution 
- + - Explicit solution 


Discharge [m3/s] 


0 50 100 150 200 250 300 350 400 450 500 550 600 
Time (Days) 


Figure 4.27 The analytical and numerical solution of the linear reservoir equation. 


Downloaded from https://iwaponline.com/ebooks/book-pdf/650788/wio9781 780400457 .pdf 
bv IWA Publishing. publications@iwap.co.uk 


Finite difference method 71 


The explicit and implicit method are applied to solve the ODE defined by (4.85) and these are compared 
with the analytical solution in order to determine how well are the numerical approximations performing. 
Application of an explicit numerical scheme to equation (4.85) gives: 


Qn = Q" 


A = -kQ" — Qn = Q" (a _ kAt) (4.86) 


Application of an implicit numerical scheme to equation (4.85) gives: 


p pu Q" HE n+l ntl — Q" 

ee ee aaa 4.87) 
Integration of equation (4.85) gives the analytical solution: 
Q0) = Q,e* (4.88) 


Table 4.1 show the values of Q(d in case equations (4.83)—(4.85) are used to determine the solution. The 
values considered for Q,, k, and At are: 


0? 


3 


Q, = 50^; k= 0017 At = 40days (4.89) 


Figure 4.27 represents graphically the results presented in Table 4.1. 


Table 4.1 Solution results for equation (4.85) for data as defined in (4.89) using 
analytical solution and explicit and implicit solution approaches. 


Time Discharge (m?/s) 
(days) Analytical Numerical solution 
solution Explicit approach Implicit approach 
0 50.00 50.00 50.00 
40 33.51 30.00 35.71 
80 22.46 18.00 25.51 
120 15.06 10.80 18.22 
160 10.09 6.48 13.01 
200 6.76 3.88 9.29 
240 4.53 2.33 6.64 
280 3.04 1.40 4.74 
320 2.04 0.84 3.38 
360 1.37 0.50 2.42 
400 0.92 0.30 1.72 
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4.6.2 PDE: Simple wave propagation 


A basic equation for the mathematical description of flood routing in a river is the conservation of mass 


equation: 
oQ oh 
ax b àr 0 (4.90) 


where Q is the canal discharge; b the canal storage width; and A the water depth in the canal. Assuming 
that the relationship between Q and A is unique (i.e., unique stage — discharge relationship)-relation (4.90): 
oQ  dQoh 


UE LUE 4.91 
Ox | dh ox Ge) 


equation (4.86) can be re-arranged as: 


dQ dh ðh 
ant oop = (4.92) 


or 


1dQ|dh | dh_ 
Rae +5 =0 (4.93) 


that simplifies to the linear advection equation: 


oh + con =0 where c= 140 (4.94) 


Using a discretization that is made up of a forward difference in time and backward difference in space, 
the solution of A(x, f) for each grid point on the x — t computational domain shown in Figure 4.28, which 
is 60 km long and has an input stage hydrograph raising from 2 m to 6 m over 6000 second. The total 
computational time is 60,000 seconds. Figure 4.29 and Table 4.2 shows the results obtained with the FTBS 
scheme for a grid size Ax = 3000 m and Aż = 3000 s; and a c= 1 m/s. Table 4.2 gives and extract for the 
first 50 minutes and for 24 km out of the 60 km long channel. 


Computational Domain| 


Boundary condition upstream 


Given B.C. 


Given I.C. 


Initial conditions 


Figure 4.28 Computational domain, initial and boundary conditions for the simple wave propagation problem. 
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Figure 4.29 Solution of the simple wave propagation. 


3000 s, over 60 km 


3000 m and At 


Table 4.2 Solution results (extract) for equation (4.94) for Ax 


000 seconds. 


long channel and during 60 


Water depth (m) 
At Distance (km) 
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(min) 
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Results are showing that in case the product c At/Ax is equal to one the advection equation 
transports the input unchanged. This product is called the CFL condition and it is described in detail 
in Chapter 6. 

In case that this product is less than 1 diffusion appears in the solution, which is due to the numerical 
approximation and should not be confused with the attenuation of a flood wave. The result that is obtained 
for a selection of Ax = 3000 m and Aż = 1500, for the same canal, time length and c = 1 m/s is presented in 
Figure 4.30. The numerical effect of diffusion can be noticed. 


Water Depth [m] 


Q = 
600 i lo 


Figure 4.30 Solution of the simple wave propagation with numerical diffusion. 


4.6.3 PDE: Diffusion equation 


The 1D diffusion equation, of a pollutant, at a river section, where the pollutant is discharged is: 


aC PC 

——— D, == 4.95 

ot m Ox ( ) 
where C is chemical concentration, and D,, is the molecular diffusion coefficient. Equation (4.95) has 
several analytical solutions depending on initial and boundary conditions. For the case of a pollutant in the 
form of an impulse input at x = 0 into a semi-infinite media (C = 0 at x = ee at all times), where C(t = 0) = 0 
at all other values of x, equation (4.95) has the following analytical solution: 


2 


C(x, t) = oe (4.96) 
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where M is the mass of chemical added at x = L, with units of mass per cross-sectional area (e.g., mg/cm’). 
This solution is valid when diffusion occurs in both directions, with the concentration profile extending 
further as time passes. 

Using a discretization that is made up of a forward difference in time and backward difference in space, 
the solution of c(x,f) for each grid point on the x — t computational domain, is given Figure 4.31, for the 


grid size Ax = 50 cm, At = 27 hours; D,, = 0.01 cm?/sec, M = 0.7 mg/cm’. The length over which diffusion 
takes place is 600 cm. 


12 
—— Initial Condition 
—*— time=140 hrs 
10 === time=300 hrs 
— -time=415 hrs 


sire time=550 hrs 


Concentration [mg/l] 


0 100 200 300 400 500 600 
Distance (cm) 


Figure 4.31 Initial condition and solution for diffusion equation. 
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Chapter 5 
Finite volume method 


5.1 GENERAL CONCEPT 


The finite volume method (FVM) is a numerical approach applicable to differential equations that represent 
various types of conservation laws. As a method it was originally developed to solve water related problems 
and has been extended to other disciplines later. FVM has the advantage that it can be used on arbitrary 
geometries and on structured or unstructured discretization of the computational domain. 

In the FDM approach, detailed in the previous chapter, the differential equation to be solved is 
written at each discretization node of the computational domain by replacing the derivatives with a finite 
difference approximation. The FDM method, however is difficult to apply when the coefficients involved 
in the differential equation are discontinuous (LeVeque, 2002). In the FVM method, discontinuities of the 
coefficients are not a problem if the computational grid is chosen in such a way that these appear at the 
boundaries of the computational domain (Durran, 1998). Finite difference approximations are sometimes 
used in FVM approach to approximate fluxes at the boundary of the control volume, the discretization 
unit of the computational domain. Hence the FVM differs from the FDM in the way finite difference 
approximation are used (i.e., for the flux of a quantity, rather than for a derivative). 

Depending on how the discretization of the computational domain is done, the FVM is also known as 
cell centred, vertex centred, box method, or balance method. 

Finite volume methods have been extensively developed by different researchers for each type of partial 
differential equation (hyperbolic, parabolic and elliptic) and dependent of the discretization type used; 
Lazarov and Makarov (1982) for Cartesian meshes, Mandel and McCormick (1989) for unstructured 
meshes; Manteuffel and White (1987) for cell-centred finite volume schemes; LeVeque (2002), Godlewski 
and Raviart (1996) for hyperbolic problems; and so on. 

Finite volume methods can be characterised as follows: 


* Problem domain is discretized into discrete elements called control volumes (CV); 

* The differential equation to be solved is integrated into a balance equations for each control volume 
of the computational domain; 

e Obtained integrals are approximated using numerical integration; 
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* The values of the unknown function and their derivatives are approximated using the values at the 
nodes of the control volume; 

* The solution of the problem over the computational domain is obtained by assembling the obtained 
discrete equations over a control volume, into an algebraic system of equations. 


An important feature of the FVM approach is the fact that the numerical flux is conserved from one 
discretization cell to another (the neighbouring cell), that is, there is local conservation of the numerical 
fluxes. This is very important especially because flux is an important element in fluid dynamics. 

The FVM principle is detailed bellow. The method is applicable for the solution of conservation laws, 
and for simplicity of the demonstration of the method the one-dimensional general form of the conservation 
law (see chapter 2) is used: 


2 fucx, DAV + [ fo «m aA - [st nav =0 (5.1) 
V A V 


Equation (5.1) shows the transport of a variable (substance, water depth, etc.) under the influence of 
advection, that is, the rate of change of variable u(x, f) inside a defined volume V together with the flux of 
u(x, f) through boundary A is the same as the rate of S(u, f). Equation (5.1) is valid at each point x and each 
time f, of the computational domain, where the conservation of u(x, f) is written. S(u, x, f) is also known as 
source term of the equation and expresses a possible volumetric exchange between conserved quantities 
u(x, f), that is, the diffusion phenomena. The term n represents the unit vector perpendicular on surface A, 
pointing outwards. 

Before showing the principle of the FVM method it is important to note that in the integral expression 
(5.1) the derivatives are zero order for the advection term (/f(u)) and first order for the diffusion term 
(S(x, £)), if it exists. The fact that the order of the derivatives is so low is important when solutions are 
sought for problems where phenomena change rapidly in space so that the spatial derivative does not 
exist, such as in the case of a hydraulic jump. Functions that are discontinuous do not have derivatives 
at the discontinuity position, hence the conservation equation expressed as a PDE, is not valid at the 
discontinuity point. Finite volume method presents the advantage that they overcome this issues, by 
considering the integral representation of the PDE and by having reduced order of derivatives, that can 
be easily treated. 

The principle of the method is illustrated on a simple example of the one-dimensional linear advection 
equation, which is obtained from (5.1) if fu) = au and S(u, f) = 0: 


du | fu) E 
(¢ TU. Jav =0 (5.2) 


V 


Consider the above equation to be solved over the one dimensional computational domain V that 
is discretized in control volumes V, defined between points x; ,, and x; (see Figure 5.1b). Volume 
integration of equation (5.2) over the control volume (CV) yields: 


Xi+1/2 Xi+1/2 


| Space J Fax = 0 (5.3) 


Xi-y2 Xi—1/2 


In Figure 5.1 points x; are equally spaced along the space domain x with a constant distance 
AX = Xian Xuan t= 1,..., N— 1. Same convention of notation holds, as in Chapter 4, that u(x;) = 
u(x; + Ax) = u;„. In Figure 5.1 control volumes are noted with V;. 
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(a) ri (b) Me 
Uj-4 T Uia Uij4 i Ui 
o ! ; | Q 
i-1 i-1/2 i Xie Xj Xv 
SSS ee 
Vi Xi V; Xi+1/2 
1D Cell centered FVM 1D Vertex centered FVM 
(c) (d) 
Í = if Nodes e e e e > Nodes 
e e E ie e 
CV : 
e cve e e 
e e e e 
2D Vertex centered FVM 2D Cell centered FVM 


Figure 5.1 Finite volume examples of control volume in one- and two-dimensional space. 


Equation (5.3) represents the integral formulation of the differential equation over control volume V; 
which is written in all CVs of the computational domain to yield a system of equations, and gives the 
solution to the problem by determining the values of the unknown variable in all CVs. 

Solution of the system formed by types of equation (5.3) requires: 


* numerical approximation of the space integrals, which can be volume and/or surface integrals, 
depending on the problem dimension; 

* interpolation, to determine the values of coefficients that describe the numerical approximations of 
the surface and volume integrals; and 

* time integration by evaluating the time derivatives. 


If u; is averaged and considered constant over CV, then the first term of the equation (5.3) can be 
integrated numerically as: 
Xi+1/2 
g” = uc Qua + X12) (5.4) 
Xi-2 
The second term of equation (5.3) integrated analytically yields: 


J 2 ax = f(u) fiu) (5.5) 


Xi-y2 


Substitution of equation (5.4) and (5.5) into equation (5.3) gives its discrete form as: 
Ui (Kis + Xiao) + fap) - fU-ar) = 0 (5.6) 


Equation (5.6) is a conservative scheme if the flux on the boundary of one cell is equal with the flux on 
the boundary of the neighboring cell. 
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The discretization of the space domain, as selected, is known as cell-centered method, where the 
unknown variables, u(x, f) are computed as average values over the CV and they are attached to the middle 
of the cell. A similar discretization is the cell-vertex method, in which the variables are averaged over a 
CV, but they are attached to the nodes of the discretization grid, that is, at the edge of the CV (see Figure 
5.1a). Figure 5.1c, d, shows two types of discretization for a two dimensional space. 


5.2 FVM APPLICATION DETAILS 


5.2.1 Step by step application of the FVM 


Based on the FVM principle illustrated above, a general step by step procedure for the application of the 
method is given for a 3-dimensional problem. 

Equations (5.1) and (5.2) are exact formulations of the problem of a phenomenon, that is, the conservation 
law principle of a quantity u — u(x, y, x, t). (Note: for a 3-dimensional function u(x, y, z, f) the notation u, 
is used when reference is given to all 3 dimensions of u). Consider the 3-dimensional advection-diffusive 
equation, which is a parabolic type of differential equation: 


du NO 9 du 
Qu (^ u^ 2d (5:7) 
i=l k 


i 


where x; are the (x, y, z) components of the space, t is time, a is the diffusion coefficient and u = u(x, y, z, f) 
is the unknown transported quantity (see Figure 5.2). 


Figure 5.2 Computational domain for the 3D advection-diffusion equation. 


The domain is divided into i (i= 1, ..., N) computational cells (Vi) — control volumes. A number of 
nodes are created at the center of the control volumes, where the unknown variable is computed. The 
balance equation for each of the defined CV is obtained by formulating the problem in an integral form, 
which is equivalent to the integration of equation (5.7). The integration over an arbitrary CV yields: 


a] PIC = T3 ndá = pos 6.8) 


i 
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where, A is the surface of the CV and n; the unit normal vector to the surface, corresponding to the x; 
component of the space. 

In the FVM approach a numerical approximation of the solution is introduced by the space integration 
of the equation, hence fluxes in space and time are to be calculated. At each CV the average value of the 
function u is known. The advection and diffusion fluxes of the equation (5.8) are calculated in three steps, 
as follows: 


(1) Determination of the values of function u for each CV: Both the value of the function u and 
its derivatives are numerically approximated for each CV. Based on these approximations both 
advective and diffusive fluxes are determined. Function u; for the CV ; (Vi) is approximated using 
a polynomial representation: 


: 
u = M a d.) (5.9) 


n-l 


where a, are multiplication coefficients of @, chosen interpolation functions (simple polynomials). 
Each 6, is valid in one of the P surrounding cells of the considered Vi control volume. Coefficients a, 
are determined from the condition that the cell averages are valid over a number of (i + m) cells, hence: 


J iav - v u 


itm™i+m? 


m=0, 1, 2,...,P—1 (5.10) 


Vitm 


where V,, ,, are the number of P cells surrounding the cell Vi. 
The unknown coefficients, a,, are determined from a system of algebraic equation, of the form: 


P 
EMO = ViemUism> with Ann = J g,dV (5.11) 


n=l V; 


i+m 


(2) Evaluation of the integrals: The integrals are approximated numerically. The most often used 
approach to evaluate the integrals is the Gauss quadrature approach. In this case the selected 
functions $, determines the number of quadrature points to match polynomials of degree P. 

(3) Temporal integration: The computed fluxes are used to advance the solution in time using an 
explicit time approximation method such as Euler, or Runge-Kutta, or the Adams-Bashforth 
class of methods. Time integration is dependent on the spatial approximation and is checked for 
stability. 


The details of the three calculation steps are further explained, along with considerations regarding 
boundary conditions and the solution of the obtained system of equations. The chapter ends with an 
example on how to apply the method for a one dimensional advection diffusion equation. 


5.2.2 Surface and volume integrals 


Consider the example of a two dimensional domain which is split in CV quadrilateral shapes. Figure 5.3 
shows the form of such a CV and the notations used on it. The nodes defining the CV are SW, SE, NE, NW 
along with the midpoints of the sides; S, E, N, and W. A special point, P, is positioned at the center of the 
CV. Figure 5.3 includes the representation of the unit normal vectors. 
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(6) yxa) 
YN 
De 
Ys 
T XE p 
l , , x(x) 
arbitrary perpendicular to the axis 


Figure 5.3 Quadrilateral control volume: (a) oblique; and (b) perpendicular on the coordinate axes. 


According to the notation in Figure 5.3, the surface integral of equation (5.8) on the left-hand-side, is 
the summation of four surface integrals, over the four sides of the CV: 


YI (« — a Zi ud (5.12) 


In equation (5.12) / is the number of sides of the CV (1.e., 4 sides), and all indices with /, represents the 
elements corresponding to side / of the CV. There are two types of fluxes; convective F^, and diffusive 
FP through the CV faces /: 


Fo = Junas, (5.13) 
Al 
and 
Ou 
F? = {fe} n,dA, (5.14) 
Al 


The area/length of side A, is the length of the CV side. For the neighboring CVs with a common side 
the total value of the flux is: 


F-F + FP (5.15) 


The sign of the flux, however is opposite for the two neighboring cells. 
The approximation of the surface integrals (5.13) and (5.14), for side / of the CV with a cell-centered 
variable, is done by steps 1 and 2, as detailed in Section 5.2.1. It yields: 


Ou (5.16) 
l 
where m; is the mass through face / of the CV. Relation (5.16) is approximated at side AE of the CV for a 


selected polynomial. AA, is the area of element /. 
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The volume integral of equation (5.9) is approximated through numerical integration. The following 
relation holds over a CV: 


V 


where Sp represents the value of S at the center of the CV (the average value of S over CV). AV represents 
the volume of the considered CV. Taking into account notations from Figure 5.3 the volume of the CV can 
be expressed as: 


1 
AV = z [eae — Xyw) One — Yaw) — Oye — Law) Oar — Yow) (5.18) 


All the integrals in equation (5.9) evaluated using (5.17) and (5.18) are added up over the whole computational 
domain. The following overall equation holds: 


Derr Ant = 25 Ney NAcy | = 


| ee (5.19) 


cv source terms 


convective fluxes diffuzive fluxes 


The values of the unknown function and its derivatives at each side of the CVs are determined using 
equation (5.19). 


5.2.3 Discretization of convective fluxes 


The convective fluxes, F^ are approximated using the values of function u at the CV center. The values of u 
at the CV center are calculated using approaches similar to FDM; upwind techniques, central differences, 
mixed techniques, and so on however they should be regarded as interpolation techniques, not finite 
differences. Interpolation is done in FVM usually in the direction of flow, that is, the assumption is made 
that the convective transport of u only takes place towards downstream, which is the case when the equation 
is of the parabolic type as is the case of the advection diffusion equation. Details on approximating the 
integrals are shown bellow on a one-dimensional equally spaced computational domain. 


5.2.3.1 Upwind approach 


The upwind method (UDS) approximates the values of u by a step function. The value of u over a CV 
on the E side of the CV (up) is determined as follows (see Figure 5.4): 
up = üp, pomme (5.20) 


Up = ug, if mg <0 (5.20b) 


A special upwind approximation is the quadratic upwind interpolation, also known as the QUICK 
method (see Figure 5.5). In the QUICK approach a polynomial is determined using an interpolation 
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through 3 points; two neighboring points P and R, and a third point, RE or PW, depending on the flow 
direction. The resulting polynomial for the approximation of up is: 


Ug = dug — djupy + (l-a + a)up, if mg >O (521a) 


ug = bup — bugg + (l1 — bj +d, )ug, if m; «0 (5.215) 


i+1 


C 
i; mg<0 


Figure 5.4 Constant approximations of u over the control volume. 


Figure 5.5 Approximation of function u using linear interpolation over two control volumes. 


with values of the coefficients given by: 


-vr a A2 Y00 -y 
BETTE ? Dey (5.22a) 
p = Vt Y, - y) p = Vele 
; EN , 2. * 
Id oe VEY Y (5.22b) 


Downloaded from https://iwaponline.com/ebooks/book-pdf/650788/wio9781 780400457 .pdf 
bv IWA Publishing. publications@iwap.co.uk 


Finite volume method 85 


and, 


y,- Xe B (5.23) 


Xp — Xp 


Relation (5.24) represents the interpolation factor for function u at side e, based on the values of u in points 
P and R (see Figure 5.6), that is: 


Lad Yelle + ad ms y Up (5.24) 


In case of a equally spaced grid: a, 2 , Ay 


8 


Figure 5.6 Interpolation at point P. 


5.2.3.2 Central Differences 


A central differences scheme (CDS) is also used for the approximation of the convective fluxes. In case of 
the CDS approach u is linearly interpolated (see Figure 5.6) and equation (5.24) applies for interpolation. 
The interpolation factor y, is defined by equation (5.23). 


5.2.4 Discretization of diffusive fluxes 


The approximation of the diffusive fluxes, FP is done by approximating the normal derivative of u at the 
CV faces. These values are expressed using the values of u in the CV centers. For example on the east face 
of the CV, Ae, a simple difference formula to express (du/dx), is: 


(à) uH Up (525) 


ox Xp — Xp 


Equation (5.25) is a central difference formula. 
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5.2.5 Evaluation of the time derivative 


Time derivatives are approximated using fractional steps (or finite differences): 


o i 
5 J udV| = V (5.26) 
V 


Vi 


5.2.6 Boundary conditions 


Boundary conditions can be specified in a number of ways; Dirichlet type, when the unknown function u 
is prescribed; Neumann type, when the flux of the unknown variable u is prescribed; and a combination 
of the previous two. Boundary conditions are equations that are added to the linear system of equations 
obtained over the computational domain in order to have the same number of equations as the number 
of unknowns. 

Consider the CV, as defined in Figure 5.3b. If u is given as a prescribed value on the w side of the 
CV, as 
U,, = u,gc then the convective flux at the boundary is approximated as: 


Fe = Muy = Hype (5.27) 


Equation (5.27) is used directly into equation (5.20). 
The diffusive flux through the boundary is calculated using the same approach as: 


(Z) zs Up — Uy = Up = UwBC (5 28) 


The algebraic system of equations obtained when the FVM approach is used has a unique solution only 
if the boundary conditions at all boundaries of the computational domain are defined. 


5.2.7 Solving algebraic system of equations 


By summing up all elements of equation (5.19) a system of equations is obtained, which gives the values 
of the unknown function at all center nodes of the CV-s. Replacing the fluxes, temporal integrals and 
boundary conditions in equation (5.19) for one of the control volumes V, (i= 1, ... , N) of the computational 
domain, the following equation is obtained: 


abuh — X,alui = bL foralli =1,...,N (5.29) 


where c represents the number of all cells around the CV, which has as a center point P. If equation (5.59) 
is written for each of the control volumes V; (i = 1, ... , N), there are N equations formed because there are 
a total of N discretization nodes. 

Schafer et al. (1993) shows how such a system of type (5.29) is obtained and solved for a one- and two- 
dimensional advection diffusion equation, in steady state (no time term). 
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Consider the steady one-dimensional advection-diffusion equation defined on domain [0, L], discretized 


into N control volumes (see Figure 5.7). Using a second-order central differencing scheme, the discrete 
equations of the phenomena, after computing all fluxes are: 


ApUp — Agly — iM, = bp 


(5.30) 
Upi! Ur uj! up 
: A i 
J | | | == | 4 4 
?— 11 E T +++ 799». 
0 PW" PÊR L 
o 9? ? 9? ? 
i] 2 ; d i M 
1 3 1 1 x 
$ i c Ll dl - 
Uy 
celli 
Figure 5.7 Discretization of the 1D computational domain. 
Values of the function over each CV, as defined in Figure 5.7, are: 
ul = ul! for all i-2,...,N (5.31a) 
u =u}! foralli—-L..,N-1 (5.31b) 
which leads to a linear system of equations, that can be represented in matrix form as: 
[a 1 Tau] bp 
ab -ah n . Up $i 
-a ab -a? 0 : d 
—ai ab -a up |=|; 
P R " bi, (5.32) 
0 -ah 
a» aP ur] 
L w xl P 
X —— | BN | 
E Vy 


u 


which can be solved to determine the values of u at each nodal point. 


Ifthetwo-dimensional case ofthe steady advection-diffusion equation is considered over a computational 
domain as the one represented in Figure 5.8, then equations (5.19) are of the form: 
ag! up! ag! ug? IBITEI ay uy ay uy — fh 


w Ow N CEP (5.33) 
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M+1 
Moe elele e 
Ui j+1 
e 
e elele e 
Uiij |Uij Uia j 
Je eus . j e e ° 
e elele e 
NN Ui j-1 
e 
1 e: e e @ |---|- e i 
0 o 
1 i N N+1 
Figure 5.8 Discretization of the 2D computational domain. 
where i= 1, ..., N are the discretization intervals on the x axis; and j= 1, ... , M, are the discretization 


intervals on the y axis. Using the notations from Figure 5.8, the matrix A of the system is: 


1,1 Ll 
üp —-aw >: 0 à —a}' 
-al? 0 
0 "m " —qN-1M 
R 
34 
A= ; ; » & A ] (5.34) 
24 
—à;. 0 
0 qua 
_AN.M N,M N,M 
| ay 0 —a, ap | 


5.3 EXAMPLE OF ADVECTION-DIFFUSION EQUATION IN 1D 


The application of the FVM approach is illustrated here on the one-dimensional advection diffusion 
equation. The domain of computation for the problem is L and the discretization of the domain is 
illustrated in Figure 5.9. In order to make the demonstration easier the numbering of CV sides along with 
the CV center nodes is done as a multiple of number 2. The CV sides are shown as ratios, while the CV 
center nodes are the numbers i — 1, ... , M (see Figure 5.9). The domain has M number of equally spaced 
volume cells. The width of a cell is Ax, which corresponds with the volume of the CV, because it is a one 
dimensional problem. 
The advection diffusion equation to be solved over the one-dimensional domain is: 


tac--a (5.35) 


Equation (5.35) has two types of fluxes; the advective flux, which is noted as F — au, and the diffusive 
fluxes noted as D = ow. These fluxes are written at any point x;, ,,. 
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Figure 5.9 Values of the unknown function over an equally spaced grid using a constant function 
approximation and a discretization of: (a) 4 cells, (b) 8 cells, (c) 16 cells. 


Equation (5.35) integrated reduces to: 
Fi-F, D,-D, 


du; By i-z iz 0$ (5.36) 
dt Ax Ax 


The evaluation of the terms F and D in equation (5.36) is done at the cell edges of the control volumes. 
Before time integration takes place the fluxes are evaluated using the step by step procedure. Values of 
function u are written as a polynomial representation, as detailed by equation (5.9): 


P —X. 
u- M as where € = E as) (5.37) 


In equation (5.37) a, are P + 1 coefficients that are determined from P + | conditions on the averages of 
the polynomials over multiple cells, as shown by equation (5.10). In equation (5.37) the coordinate system 
is centered to the cell and scaled so that the cell length is from —1/2 to 1/2. 
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Applying equation (5.11) coefficients are computed as: 


x n-l 


Ann = J EE dx forn=1,2,...,P+1 (5.38) 
X | 
X qo— X3 i X 47X | i 
B Ax; me. i3 Hy i-z 
== Ax Am (5.39) 
1751 x (5.40) 
p " 2 di s; — "mti z 
SAN n y um Sl ~ Ax; 


Values of A,,,, are determined and the system of equation (5.38) is solved. Different polynomial 
representations for the function u are considered and shown bellow, for an initial condition of a cosine 
function. 


5.3.1 Constant unknown function 


The constant value of u, as shown in Figure 5.4, is the simplest case, in which the approximation is given 
by the relation: 


u = dg (5.41) 
In such a case a, is the only unknown coefficient, hence only one constraint has to be fulfilled: 


X 1 
ity 


adx = uj^x; 
: (5.42) 


Equation (5.40) states that the volume integral over the control volume is equal with the function value 
times the control volume length, hence: 
d = uj (5.43) 


Two approximations are available for (5.41), either the left value of u, in cell i, or the right value of u in 
cell i+ 1. The flux is advective, hence the information usually is traveling from upstream to downstream, 
therefore, the solution is: 


2 
u= 5.44 
it Un fa, <O SER 

m 

2 


The diffusive term cannot be computed when constant functions are used because the derivative of a 
constant is Zero. 
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5.3.2 Linear variation approximation of the unknown function 
A better approximation than the constant one is the linear variations of the unknown function over two 
adjacent control volume cells. If a linear approximation is considered for u then: 


u = a + ač (5.45) 


In equation (5.45) there are two unknown coefficients, ag and a,, hence two conditions are set in order to 
lu aegra them. Consider two consecutive control volumes ; and (i + 1) so that X SxS, m ; hence 
-4<&< 3. The two necessary conditions to find ay and a, s are: 5 


1 


J (a + ag)dg = u; (5.46) 
D "n 
fea + add = eu, (5.47) 


After performing the integration, the following system of equations is obtained: 


1 Oa) “i 
For an equally spaced discretization, it yields: 
Ay = Uji, A, = UW — Uj (5.49) 
The linear interpolation of u over the two adjacent cells gives: 
Uu = U; + (ua — Uj) (5.50) 
hence 
u=(1-€)u, + uzi, for- j «és ; (5.51) 


The unknown function u at the edge of the control volume is obtained for values of € = +5. 


5.3.3 Parabolic variation approximation of the unknown function 


Considering an equally spaced discretization grid and a parabolic interpolation of the unknown function 


u, yields 
_ | Ua + 26u; — uj , Uis — Ui- uj, — 2u; + uj 2 
=| 24 ll 2 lel ? é (5.52) 
The interpolation is centered on the control volume cells (i — 1), i and (i + 1), and as in previous cases: 
eee" £22 4 (5.53) 
2 2 Ax 
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At the edge of the control volume, where x;,,,. and € = 1/2 the value of the unknown function u is: 


—u; 4 + Su; + 2uj,, 
u= 


| : (5.54) 
its 


In order to assess the how good are the different forms of representation of the unknown function an 
example is considered, for which the analytical solution is known. The difference between the two results, 
the analytical and numerical solution, is the error, which depends on the discretization. The more refined 
discretization, the smaller the error. 


5.3.4 Error of the approximation 


Let u be coszx over the interval —1 € x < 1, as the initial value. The advection equation solution at all time 
steps is the same as the initial value as it was demonstrated in Chapter 2. The advection equation transfers 
the initial profile unchanged along x axis. 

The result of the approximation using constant, linear and parabolic representations of the function u 
are shown in Figures 5.9, 5.10 and 5.11, respectively, and the error at each cell edge is given by the height 
difference between the solid curve (analytical solution) and the square symbol. 


-1 -0.75 -0.5 -0.25 0 025 05 0.75 1 -1 -0.75 -0.5 -0.25 0 0.25 05 0.75 1 


-1 -0.75 -0.5 -0.25 0 025 05 075 1 


Figure 5.10 Values of the unknown function over an equally spaced grid using a linear function 
approximation and a discretization of: (a) 4 cells, (b) 8 cells, (c) 16 cells. 
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Figure 5.11 Values of the unknown function over an equally spaced grid using a parabolic approximation 
and a discretization of: (a) 4 cells, (b) 8 cells, (c) 16 cells. 


In Figure 5.9 the representation of the solution using constant approximation and a number of 4, 8, and 
16 discretization cells is given. The solid curved line is the exact solution, while the dotted lines represent 
the approximate solution. Similarly, Figures 5.10 and 5.11 shows the representation of the solution for the 
case of linear and parabolic approximation, respectively. The representation of the linear approximation is 
done over two cells, while the parabolic one is represented over three cells. 

Analysing the results it can be noted that cosine wave is not well represented in case that a coarse 
discretization grid is used along with a constant approximation for the unknown function. In case that a 
linear and/or parabolic approximation is used with a finer discretization grid, then the solution is better, 
however there are still several values that are out of bound. 
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Chapter 6 
Properties of numerical methods 


6.1 PROPERTIES OF NUMERICAL METHODS 


Each solution method for an ODE or a PDE using numerical approaches gives an approximate value for 
the unknown function. The errors in the solution come from the way the derivatives are expressed on the 
discrete parts of the domain. The main question that arises is how good the approximation is? how close 
is the computed solution to the analytical solution of the differential equation. This can be mathematically 
achieved by imposing the condition that the numerical solution converges to the exact analytical solution, 
however this is not an easy task because the analytical solution is in general not possible to be determined. 
The study of the convergence is actually done by looking at the error behaviour. The study is based on 
the theory of truncation error (the error between the analytical and approximate solution), which states 
that refinement of the discretization grid gives better numerical solution (more accurate) (Abbott & Basco, 
1989; Fletcher, 1991; Guinot, 2003). 

In practice an indirect procedure is used, by imposing that the system of algebraic equations obtained 
through numerical solution approach, to be consistent with the governing differential equation. In addition 
the numerical solution should be stable. Consistency and stability ensures convergence. This is stated by 
Lax theorem of convergence, which is presented bellow. 

The process of discretization and ensuring convergence of numerical solution to the analytical solution, 
as defined schematically by Fletcher (1991), is presented in Figure 6.1. 


6.1.1 Convergence 


A solution obtained by a numerical scheme is by definition convergent if it is closer and closer to the 
analytical (exact) solution of the differential equation (ODE or PDE) when the grid size of the discretization 
tends to zero. The difference between the exact solution and the value of the solution obtained with 
numerical approximation is called solution error. 

There are two different meanings of convergence in numerical analysis; if the discretized interval is 
getting finer and finer after discretizing the continuous problems, the solution is convergent to the true 
solution; and for an iterative scheme, convergence means the iteration gets closer and closer to the true 
solution as computation progresses. 
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PDE describing a 
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Consistency + Stability ————————- Convergence 


Figure 6.1 Relationship between convergence, stability and consistency. 


In order for an approximate solution to converge to the exact solution of a discretised equation (ODE or 
PDE) as the time and/or space step sizes goes to zero, two conditions must be met: consistency and stability. 


6.1.2 Consistency 


A solution obtained with a numerical scheme is consistent if it gives a correct approximation of the 
differential equation as the computational, of time and/or space, step is decreased, which means that local 
truncation error goes to zero as step size go to zero. Consistency is usually verified using Taylor Series 
expansion. The analysis of the consistency of a scheme is done for each particular studied scheme, and it is 
further presented in this chapter for the schemes that were introduced in Chapters 4 and 5. 


6.1.3 Stability 


Stability refers to the magnitude of the change of the error. A scheme is stable if any initially finite 
perturbation remains bounded as time grows. A scheme is stable if initial errors or small errors at any time 
remain small while computation progresses. It is unstable if initial errors or small errors at any time get 
larger and larger, or eventually get unbounded, as computation advances in time. 

Though it may seem a very simple issue, it is interesting to mention that stability is an important 
property of the numerical methods which is not easy to notice and has not even been discovered by the 
one who defined the finite difference method for the first time, L.F. Richardsond in his book on weather 
predictions using numerical methods, which was issued in 1922. The issue of stability was noticed by 
Courant, Friedrichs and Levy in late twenties (hence the name CFL condition); was well defined by von 
Neumann in 1940’s and rigorously defined by Lax and Rychtmyer in 1950’s. 

Numerical methods can be classified in absolutely stable, conditionally stable and unstable. 

Stability verification is more complex than consistency. Several methods are available to carry out 
stability analysis: 
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* Matrix method, which is based on finding the eigen values of the matrix representation of the solution 
given by the numerical scheme. It is used mainly for FDM method. 

* Fourier method, or also known as von Neumann analysis, in which the solution given by the 
numerical method is expressed as a complex Fourier exponential. The substitution of the Fourier 
representation of the solution into the original equation allows for the analysis of growth or decay of 
the error. The method however is limited to linear equations. 

* Domain of dependence, in which the domains of dependence of the analysed differential equation 
and the numerical representation of the solution are compared. 


6.1.4 Lax's theorem of equivalence 


Forthe restricted class of linear problems Lax's theorem states the relationship between consistency, stability 
and convergence for finite difference approaches, however the theorem is applicable for any numerical 
approach that involves discretization of the computational domain and representation of the differential 
equation on such a domain (Richtmyer & Morton, 1967). The theorem states that consistency and stability 
are a necessary and sufficient condition to obtain convergence. Any of the two conditions alone are not 
sufficient to have convergence. Fletcher (1991) states that for the case of fluid problems the theorem cannot 
always be applied rigorously and should be looked at with care, because fluid problems are non-linear in 
nature and have mixed boundary and initial conditions. 

Lax theorem is a consequence of the analysis of the error between the numerical solution and the exact 
solution of an equation at a node of the Discretization grid. The analysis of the error, which is expressed as 
the difference between the exact and the approximate solution, shows that the error at time level n depends 
on the error introduced at initial time level; on the truncation error; and on the discretization grid. The 
approximate solution can be stable if the error is bounded and the truncation error is made very small, 
which is equivalent with the consistency condition. 

In practice the above analysis on the bounding of the error is expressed by putting the condition that 
the growth of the error of the numerical solution should be bounded at all times, which reduces to the 
so-called Von Neumann stability criterion. The criterion is expressed as: 


lA < 1 (6.1) 
where A is called amplification factor and it has the expression: 


yt! 


A= (6.2) 


n 


u 


In expression (6.2) u” is the numerical approximate solution obtained for the unknown function, at time 
level n of computation, and u”*! at time level (n + 1). 


6.2 CONVERGENCE OF FDM SCHEMES 
6.2.1 Convergence for ODEs 
6.2.1.1 Consistency 


In order to analyse consistency of a numerical approximation the truncation error is analysed. In case of 
ODEs, as defined by equation (2.21), an analysis of the error obtained from time level n of computation to 
time level (n + 1), gives: 


TE (At) = u(t” + At) — u™! (6.3) 
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in which uz"! is the approximate solution computed at time level ( + 1) and u(t" + Af) is the exact value of 
the function u(f) at the node of the grid where the evaluation of the error is carried out. If u(f) is placed in 
equation (6.3), the error is a function of only one variable, Ar. Using Taylor series expansion to express the 
value of TE(At) around zero, yields: 


are] ER (At) m 
0 0 0 


TE (At) = TE (0) + At d(At) 2! d(Aty 3! d(Atp 


(6.4) 


According to equation (6.4) any numerical method is a good approximation if the value of the error is 
as small as possible or equal to zero, hence the smaller the Ar the smaller the error. Equation (6.4) shows 
the magnitude of the error of the approximation. 

A numerical approximation method is defined to be of order p, if the first p terms of the equation (6.4) 
are zero. Hence a method of order p generates an error for which the first term is: 


TE(At) = 


pti pti 
(At) d mp] 65 
0 


(p D! d(AryP*! 
Consider the weighted method for approximating the solution of an ordinary differential equation. 
According to the weighted method the approximate solution, as shown in Chapter 4, is: 
u™*! = yu" + At-(1—6)- f(u",t") + At- 0- f(u, t) (6.6) 


Equation (6.6) is also known as the multi-step method solution. Ascher and Petzold (1998) defined a 
general form for this solution as: 


k k 
Nus a = At: 2p, à dna (6.7) 
m=0 m=0 
where: 


— Qt, D, — coefficients of the method; 

— u” — values of the function u(t) at grid point n (see Figure 4.7); 
— f” — value of the derivative of u at grid point n; 

— k — the number of past steps applied to solve the ODE. 


For the particular case of equation (6.6) the values of the coefficients in equation (6.7) are: 
a =1 o, --L B=0 Bj21-06; k=l; (6.8) 
The error of approximation for (6.7) at time level n is: 
k k 
TE(At)], = by ag — At- $5 fem (6.9) 
m=0 m=0 
The order of the method can be determined by analyzing the error, which according to (6.4) is: 


n P 
TE(At) = VC, AD- a (6.10) 


p=0 
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where the coefficients C, have the expression: 


C, = Yo. 


m=0 


(6.11) 
E 1 č 1 k 
C. = (-D^-—- Xm -Qn t: m" . p, 
ý 2 p! 2, (p — D)! È p 
Calculation of the above coefficients determines the order of the method, that is, if Cj, C, ... Ton =0 
and C,,, # 0 then the method of approximation is of order p. In case of equation (6.6) these coefficients 
are: 
Cc, =1-1=0 
C, 21-0-(1-0)20 
1 
C, =-5+(1-6)= 05-6 (6.12) 
1 1 1 80 
Ge) 8) t, 


Based on (6.12) both explicit and implicit methods are first order. 


6.2.1.2 Stability 


Stability means that if an initial small error is introduced in the initial conditions of the ODE, this should 
remain bounded, hence the computed solution at the next time level should only slightly change with 
respect to the original solution. Unstable methods produce oscillatory solutions. If the growth of numerical 
error is not bounded then the solution is unstable. 

Stability condition for an ODE is usually checked using Fourier analysis and the condition expressed by 
equation (6.2), or what is known as Von Neumann method of analysis. 

Ascher and Petzold (1998) show that multi-step methods can be defined using two characteristic 


polynomials: 
k 
pZ)= Va, ze» (6.13a) 
m=0 
k 
o(Z) = pas ‘ VAR (6.135) 
m=0 


Characteristic polynomials (6.13) are expressed in Z, which is a complex number. For equation (6.6) 
these two characteristic polynomials have the form: 


p(Z)=Z-1 (6.14a) 
o(Z)=0-Z+(1-8) (6.14b) 
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The numerical solution is stable if the ratio of these two polynomials is less than 1 (Szymkiewicz, 2010), 
which is equivalent with condition (6.2), hence: 


Pom pO) pe) 


EZ) oP) (6.15) 


where i is the imaginary unit, and @ is the argument (angle) of the complex number, which takes values in 
the interval (0, 27). Knowing the trigonometric expression, Euler, for a complex number: 


e? = cosQ t i- sing (6.16) 


condition (6.15) is: 


B Z-1 B e? —] _ cosy +i-sing —1 
 6-Z+(1-0) 6-e% +(1-6) O(cosp +i- sing) + (1— 0) 
(cosg — 1) + i- sing (6.17) 


7 (0 -coso -1—0) t i- OsinQ 
A complex number has a real and an imaginary part and A as a ratio of two complex numbers is: 
A = Re(0,9) + i - Im(0,Q) (6.18) 
with 


(cos — 1(0 - cosQ + (1+ 0)) + 0 - sin?(g) 
(0 - coso + (1— 0) + €? - sin? (Q) 
(0 - coso + (1— 8))- sing — 0 - (cos@ — Ising 
(0 - coso + (1 — 0)? + 0? - sin? (o) 


Re(0,9) = 
(6.19) 
Im(0,9) = 


Based on relations (6.19) the region of absolute stability can be represented graphically in the imaginary 
plane (see Figure 6.2) for different values of 0. Weighting factor 0 gives the extent of the region of absolute 
stability, which is the region for which equation (6.1) holds. 

In case of the explicit method 0 = 0 and region of stability is very small, while the maximum extent for 
stability is obtained for 0— 1, which is the case of the implicit method. There is one family of numerical 
methods for ODEs that are called A-stable methods. By definition (Szymkiewicz, 1993, 1995, 2010) a 
difference method is considered A-stable if its region of absolute stability covers the entire left pane of the 
complex plane. 

Another important term related to stability analysis of an ODE is the phase error. The phase error is 
the difference in phase between the analytical and numerical solution and is determined by analysing the 
argument of the amplification factor. The amplification factor defined by equation (6.2) is expressed in 
complex form as: 


nel 


v= IAl- ein (6.20) 


u n o 


A = 
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Figure 6.2 Region of stability for multistep methods of solution for ODEs. 


In equation (6.18) the term u is the argument of the amplification factor. If the analytical change of 
phase for an ODE is noted as u,, the ratio of the numerical and analytical phase is called the relative phase 
error and is expressed as: 


R=% (6.21) 


Value R> 1 means that numerical phase exceeds the analytical phase, hence scheme is called 
accelerating; while R < 1 is called decelerating scheme. Further detailed analysis for accuracy and stability 
of ODEs can be found in Dahlquist & Bjorck (2007) or LeVeque (2007). 


6.2.2 Convergence for PDEs 


As it was the case for the ODEs, a numerical solution of a PDE is a good solution if it is convergent to the 
analytical solution. The approximate solution tends to the exact one, if the grid dimensions are reduced. 
In case of a one dimensional problem, in the (x, ?) computational space this is equivalent with reducing 
systematically the dimensions of At and Ax (i.e., u — u(x;,t"), as Ax, At > 0). 
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Convergence for PDEs is studied, similarly as for ODEs, by looking at the difference between the exact 
and the numerical solution of the partial differential equation, at the same discretization grid point. In case 
of a one dimensional problem the error is expressed as: 


TE? = u(x,,t") —u? (6.22) 
in which u(x; t") is exact value of the function u(x, f) at the grid point (x; t"), and u? is the approximate 
value of the function u(x, f) at the same grid point. 

In general it is very difficult to proof that the numerical solution converges to the exact solution of 
a PDE, even for simple cases (Fletcher, 1991). The Lax theorem of equivalence is also used for PDEs, 
though it is applicable just for a restricted class of problems, the ones that are linear. All fluid problems 
are in general of non-linear nature, hence Lax theorem is used as a necessary condition for convergence 
and as a useful indication for excluding inconsistent approaches. Each problem has to be treated on a case 
by case basis. 

For example in case of diffusion equation an analytical solution is available, hence a numerical 
approximation of the solution can be obtained on a successively refined grid. The error of approximation is 
computed for each approximation. If the error reduces to zero, as the grid discretization is reduced implies 
that approximate solution is convergent to the analytical solution. 

The exemplification of the procedure to demonstrate convergence of a numerical solution to 
the analytical one is shown for the case of the linear advection equation, as represented by equation 
(4.43) in Chapter 4. A one dimensional linear advection equation with positive constant advection 
velocity a is: 


DE utto (6.23) 


In case of an equally spaced grid, on both x and : directions, the forward time backward space (FTBS) 
up-wind numerical solution to the equation is given by equation (4.48), as: 


utt! = Cr-ut,+(-—Cr)-u? (6.24) 
Recall that Cr is the Courant number: 


At 
SS 6.25 
Cr=a ( ) 


Consider that u(x, f) represents the water depth in a 10 km long canal, which has an initial water level 
of ug = 5 m. Based on the canal geometry the advection velocity is a= 0.5 m/s. 
For the given initial condition of the form: 


u O<st<tr" 
0 tm? 


t 
w-(2-5), t" <t<2t” 
i" 


0, t 2 21" 


(6.26) 
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equation (6.23) has an exact solution as shown in Figure 6.3. The exact solution at a later moment in 
time and space is the same as the initial condition, that is, initial condition is, advected unchanged 
towards downstream. For different values of grid discretization, such that Courrant number is kept 
constant, Cr — 0.5, the result is presented in Figure 6.3. Comparing numerical solutions with the 
analytical solution, it is noticed that the numerical solution tends to the analytical solution as the grid 


dimensions are reduced, which shows that the up-wind FTBS scheme converges to the analytical 
solution. 


6- 
95r x-0 x-2500m 
5| Boundary Condition — — ^, Cr=1, At=100 sec 


Cr=0.8, At=80 sec 


Water Level (m) 


__....Cr=0.5, At=50 sec 


‘| Cr=0.5, At=100 sec 


Pe 


000 2500 3000 3500 4000 4500 5000 
Time (sec) 


0 500 1000 1500 2 


Figure 6.3 Linear advection equation's initial condition, exact solution and numerical approximations by 
the FTBS up-wind scheme. 


The above demonstration shows that for the particular case of advection equation a particular numerical 
scheme converges to the analytical solution, however this is not easy to do for other types of equations. 


6.2.2.1 Consistency for PDE 


A numerical approximation of a PDE is consistent with the equation itself if in the limit, when the 
discretization's dimensions tend to zero, the approximation becomes the PDE at every node of the grid. 
Consistency is tested by replacing the values of the unknown function by their Taylor series expansion 
around the considered node. The equation obtained by using Taylor series expansions is not the same 
as the original PDE, it contains the original equation plus additional terms (reminder). The condition 


of consistency requires that the structure of the reminder is as such that it tends to zero as the grid is 
refined. 
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In case of the numerical scheme (6.24) for the linear advection equation Taylor series expansion of the 
point values u/ ,, and u around point (i, n) are: 


n 


ou Ar Qul — Ax? Bu 
uly — u; Ax a + 2 ag 6 pu oe (6.27) 
DT ntl AP An n+l AC PET n+l 

nn (6.28) 
u; d AP y 2 an, 6 a| 
Substitution in equation (6.24) gives: 

ou ou 

a IR (6.29) 
X +a E» TE(Ax, At) 

where TE(Ax, At) is the reminder term (or the truncation error term): 
At u  aAx u At Bui as Hu 
= : 3 . ] - (6.30a) 
Ro 3 QE. a € UB 
Rearranging (6.30a) gives: 
2 
a ndis dod) aes 2 Z O(Ax2) + OA?) 
ot ox 2 g? —— —— (6.30b) 
TE2 TE3 
TEl 
TE(Ax, At) 


where O(Ax?) and O(A?’) are terms that contain higher order derivatives starting with third order derivatives 
with respect to x and t. 

In equations (6.29) and (6.30) the indices i and n are omitted, however expressions are valid for the 
discretization node (i, n). The term TE(Ax, Ar) is dependent on the selected approximation method. As the 
dimensions of the grid discretization tend to zero TE(Ax, Ar) tends to zero and consequently equation (6.29) 
becomes the original PDE (equation (6.23)). 

The approach demonstrated to prove consistency for the linear advection equation can be applied for 
any numerical scheme. 


6.2.2.2 Stability for PDE 


A numerical solution for a PDE should not only be consistent with the discretized PDE, it should also be 
stable. Stability is the property of a numerical method that keeps any initial error bounded. Errors may be 
due to boundary conditions, parameters, or due to round-off errors of computers. A numerical method is 
considered stable if a small variation in the input data causes small variation in the final result. If a method 
is not capable to keep initial errors bounded and they are growing in time and space then the method is 
called unstable. 

Considering again the example of the advection equation and further approximating the solution, using 
values for At and Ax such that Cr » 1 (see Figure 6.4), it is noticed that solution is affected by oscillations 
and does not represent accurately the exact solution. Oscillations increase in amplitude over time, as well 
as in space. 
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Figure 6.4 Linear advection equation's numerical approximations by the FTBS up-wind scheme for 
Courant numbers greater than 1. 


Conditions of stability of the advection equation are investigated bellow. According to equation (6.22) 
each approximated value of the unknown function at a discretization node is a difference between the 
exact value and the error. For the three terms in equation (6.24) it yields: 


up E u(x,,t”*!) = TE 


(6.31) 
u; = u(x;,t") - TE? 
ui = u(x,t”) — TES, 
Relations (6.31) are replaced into equation (6.24) to obtain: 
u(x;, t1) - TE"! = Cr(u(x,_,, t") — TE?) + (0 — Cr)(u(x,, t") — TE?) (6.32) 


The exact solution of the linear advection equation satisfies relation (6.24) of the numerical scheme: 


u(x;,t"*) = Cr- u(x,t”) + (0 — Cr): ux; t”) (6.33) 
hence: 
TE?" = Cr- TE}, + (1 - Cr): TE? (6.34) 
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Equation (6.34) is written for all the points of the discretization grid. The analysis using von Neumann 
method implies to expand all errors from time level n as finite complex Fourier series and to analyze how 
a single component of the series behaves from time level n to time level (n 4 1). If each component of the 
Fourier series is stable then the numerical solution is stable. If one single component of the Fourier series 
is not stable then the solution is not stable. 

Based on the above statements looking at the error term of time level n, TE? expressed in terms of 
Fourier series (Abbott & Basco, 1989) gives: 


KF 
TE? = X ao ; At) , gi Fm jx (6.35) 


k=1 


(Note: Throughout the book indices notation i is used for space discretization. In order to avoid confusion 
with the unit imaginary number, j is used as an indices notation along x axis, in the Fourier series (6.35).) 
In relation (6.35) notations are: 


— j: 2,3,..., M, — index of the discretization nodes, 

i: the imaginary unit, 

— n: index of time level, 

k : order of the Fourier component, 

— m: wave number, 

— Ax: space discretization interval, 

Ao : Aj: Fourier amplitude of component kth at the time level n, in which Ao is the initial amplitude 
and A; is the amplification factor from one time level to another, 

— KF:the finite number of considered Fourier components. 


Equation (6.35) corresponds to a sine wave travelling along x axis with an amplitude (Ao) multiplied by a 
term (A7) at each time step. Such a wave is called a harmonic component, or Fourier component, or Fourier 
mode. For the case of linear PDE with constant coefficients, the Fourier modes behave independently one 
from another, being amplified or damped based on their wave length only. 

The relation between wave number (m) and wave length (A) is: 


E 


m=- (6.36) 


If parameter N is introduced as being the number of discretization grid intervals available in one wave 
length: 


A ED. (6.37) 


and equation (6.36) is multiplied by Ax on both sides then the term @ (the phase of the error) called the 
wave number, is defined as: 


mys Ape 6.38 
go =m- Ax N (6.38) 


The shortest waves represented at the considered grid point is 2Ax, while the longest one is ee, hence 
2 < N < œ, which gives the range of the wave number in the interval [0, 7]. 
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The error term (6.35), for only one Fourier component, k = 1, is: 
TE; = A". eii (6.39) 


Substitution of equation (6.39) into (6.34) yields: 


A™! . ef @i = Cr. An | g?U-D + (1 — Cr)A" - eh? (6.40) 
The ratio: 

Annu . 
A= > 1-Cr+Cr-e'? (6.41) 


is the amplification factor of the k-th component of the Fourier series representation of the error function 
TE while passing from time level n to time level (n + 1). The error does not increase if the modulus of the 
amplification factor is less than one, that is, IAI < 1. 

Knowing from (6.16) that e^? = cos o — i - sing, yields: 


A-l-Cr-«Cr:cosQ- i- Cr. sing (6.42) 


Equation (6.42) shows that A is a complex number, hence condition on the modulus is: 


IAI = J0 — Cr(1 - cos)? + (Cr - sing? <1 (6.43) 
Considering the trigonometric relationships: 


9 9 9 oP 


= cos? + — sin? — ing = 2sin = cos — (6.44) 
coso = cost ^ — sint 7 and sing = 2sin 7 COS 
Equation (6.43) gives: 
IAP = 1- 4Cr(1 - Cr)sin? 7 (6.45) 


The wave number q ranges in the interval [0, 7]. For @=0 the inequality (6.43) is verified (1 = 1). In 


case @= z the maximum value of sin? ($) in expression (6.44) is 1. The term Cr(1 — Cr) is a parabola 


which has a maximum value of 0.25 for Cr = 0.5. Hence the amplification factor remains less or equal to 
] as long as: 


0<Cr<l1 or ws (6.46) 


Equation (6.46) is known as the Courant-Friederichs-Levy condition (CFL condition) and is valid for 
hyperbolic equations solved by explicit schemes. The CFL condition is a very important one, however 
is not a sufficient condition for convergence. It states that a necessary condition for convergence of a 
numerical scheme to the analytical solution is that the numerical domain of dependence contains the 
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analytical domain of dependence. Courant number is also known to measure the ‘numerical speed’ of the 
numerical solution. 

The FTBS scheme for the solution of the linear advection equation proofs to be a consistent and a stable 
scheme, if Cr « 1, hence it is a convergent scheme. 


6.2.2.3 Numerical diffusion and dispersion 


Consistency and stability addressed the issue of convergence, however they are not the only measures 
for the good fit of a numerical solution with the analytical one. There are other errors that can occur 
while solving a PDE, which can be of equal importance. For hyperbolic type of equations, such as the 
case of the advection equation, the solution has the form of a wave, hence wave dynamics are important. 
It is important that the solutions of such equations are not distorting the wave propagation. The applied 
numerical solution should not change the wave amplitude nor its phase celerity, which is equivalent with 
saying that the numerical solution should not generate any dissipation or dispersion of artificial nature. A 
numerical method, which produces artificial dissipation is called dissipative and a method which generates 
artificial dispersion is called dispersive. 

Numerical dispersion and dissipation of the solution of hyperbolic and parabolic PDEs are more 
difficult to overcome as compared to elliptic equations. This is due to the fact that the equations contain 
advective terms. The effect of numerical diffusion can be noticed on Figure 6.5, which shows the effect 
of varying numerical parameters on the solution of the advection equation (6.23). In Figure 6.5 a one- 
dimensional Gaussian shaped initial condition to the advection equation for which a= 1 is shown. 
Numerically the same solution as the analytical one is found for Cr = 1, but for Cr = 0.5 the initial shape 
is damped. Numerical diffusion causes amplitude error, and tends to make computed profiles smoother 
than they should be. In case of transport of the pollutant in a river numerical diffusion leads to an 
underestimation of peak values in the concentration profile, which is undesirable. Diffusion problems 
are related to the stability of a scheme and are produced by the first term in the expression of TE. 


1.2r 
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D : : i 1 ^s 
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Figure 6.5 Numerical diffusion. 
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Because numerical diffusion is not a desirable phenomenon, a better scheme to discretise the advection 
equation is usually used. For example the Preissmann scheme with 0 = y= 1/2 is second-order accurate, 
which implies that the first term in the truncation error is zero. The equation is consistent with the 
following one: 


Qu, du ðu 52 


Equation (6.47) is in fact a dispersion equation. Instead of smoothing profiles artificially, a dispersive 
scheme tend to make them sharper and create oscillations. This is the numerical dispersion effect. Figure 
6.6 shows a typical example of computed profile using a dispersive scheme. 


1.2, 
t=80 sec 
t=0 i Cr=0.5 
A Cri 
0.9 -+ 
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[e] 
[e>] 


0.3 F 


0 10 20 X30 4D 60 60 70 80 90 100 


Figure 6.6 Linear advection equation's numerical approximations by the a dispersive scheme. 


A typical example of the effects of numerical dispersion is that the computed concentrations of a 
pollutant in a river may become negative. Dispersion is due to errors introduced when estimating the 
derivatives. Numerical dispersion causes phase errors. 


6.2.3 Amplitude and phase errors 


As it can be seen from equation (6.42) the amplification factor (A) for the advection equation is a complex 
number. Both its amplitude (IAI) and its argument (u) characterise the numerical scheme (Note of caution: 
@ is the phase of the error related to wave number, not the argument of the amplification factor). The 
modulus of A shows how the various Fourier modes are being amplified from one time step to the next, 
while its argument (or phase u) show how fast the modes are travelling. From equation (6.45) it can be 
noticed that various Fourier modes are not being amplified in the same way from one time step to the 
next, because for different wave lengths correspond different values of mAx. Similarly Fourier modes do 
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not necessarily travel at the same speed. The graph that shows how the Fourier components are amplified 
is called amplitude portrait and the one showing at what speed the Fourier components travel is called a 


phase portrait. 
For the FTBS numerical scheme for the advection equation the argument of the amplitude factor given 


by (6.42) is: 


Im —Cr - sing 
= = aie S 6.48 
U = arg(A) actan( ZD) arctan( ET ess] (6.48a) 


Following the same approach as for the FTBS scheme the argument of the exact solution results in: 


U, = arctan(-Cr : 9) = -Cr - m: Ax = c d (6.48b) 


Based on (6.48) the ratio of the numerical to analytical phase speed R (see equation (6.21)) is calculated. 
Amplitude portrait is shown in Figure 6.7 for the advection equation for the FTBS scheme for different 


values of the Courant number. 


12r 


0.9- 


0.6- 


Amplitude factor 


0.3- 


Figure 6.7 Amplitude portrait for advection equation using FTBS numerical approximation. 


Figure 6.8 shows the decrease in amplitude of Fourier components for the FTBS scheme after one time 
step and after 100 time steps, using a value of Courant number of 0.5. Analysing the amplitude decrease, 
it can be concluded that FTBS scheme has strong numerical diffusion. This is due to the fact the scheme 
is of first order and truncation errors are large. FTBS scheme is not recommended for the solution of the 
advection equation, except for cases when boundary conditions are outflows. 
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Figure 6.8 Amplification factor of the wave components for advection equation. 


6.3 CONVERGENCE OF FVM SCHEMES 


6.3.1 Convective fluxes 
Convergence of the FVM schemes are analysed for the discretization of convective and diffusive fluxes, as 


they have been developed in Chapter 5. 
The discretization of the convective fluxes in equation (5.19) was done using the upwind and the central 


schemes. 
In case that an upwind scheme is used the Taylor series expansion of the solution of u, at up, around the 


middle point P (see Figure 5.7) gives: 


ou (xg — Xp) EJ 
Up = up ct (Xp +x ( + +TE (6.49) 
E P E P ox 2 2 ox s 


In equation (6.49) TE is referred as previously as truncation error. From equation (6.49) it can be seen 
that the method does not dependent on the computational grid. The interpolation error of the scheme is of 
Ist order. The first term in the TE of the convective flux F£ is: 


TE, = m,(x, — Zt (6.50) 
P 
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The error as it is expressed by equation (6.50) is called numerical diffusion (sometimes artificial 
diffusion), because the term is like a diffusive flux. The advantage of using an upwind scheme is that such 
a method generates an unconditionally bounded error. The magnitude of the numerical diffusion is given 
by the coefficient m,(x, — x,). In case that the transport takes place perpendicular to the control volume 
face, then the approximation of the convective fluxes is good, because the term (du/dx), is small in such a 
case. In case that the condition of transport being perpendicular on the face of the control volume is not 
met then the approximation may be inaccurate and for large values of fluxes (such as large velocities) it is 
required to use very fine grids in order to obtain good accuracy of the solution. Very fine grids would make 
the difference x, — x, very small. 

In case that a central scheme is used for the discretization of convective fluxes then the approximation 
has an interpolation error of a second order, as it is illustrated by the interpolation in equation (5.12). The 
second order of the interpolation error is valid for both an equally spaced grid and an unequally spaced 
grid. Taylor series expansion of u, around the middle point P gives: 


"NT (à egeo 2 (6.51) 
P P 


Evaluating expression (6.51) at x, and x; and comparing them leads to: 


= ET 2 
U, = y,ug * (01 y up (x, mex 2E tT, (6.52) 
P 


which shows that the first term in the error depends quadratically on the grid discretization. 


6.3.2 Diffusive fluxes 


The analysis of the error magnitude of the diffusive fluxes is done by analysing the difference of the Taylor 
series expansion at points x, and x,, around point x,. Hence: 


(à) _ Me = Up cg mcr [ex [d ATE (6.53) 


ox Xp — Xp 2(Xg — xp) ox? 6(Xxg — Xp) ox? 


Analysing expression (6.53) it can be noticed that for an equally spaced grid the second term of the 
expression becomes zero, hence a second order error results. In case that the grid is not equally spaced, 
shows that the term is proportional with the expansion rate of the neighboring grid cell, €: 


(= &Xx, - 263 | with €, =T% 


(6.54) 


2 X, — Xp 


which is equivalent with stating that the larger the expansion, the larger the error. This gives modellers the 
indication on how should unequally spaced grids be build, and that expansion from one grid cell size to 
another should be done smoothly and with care. 

There are other possibilities available for approximating derivatives when using finite volumes, such 
as forward or backward differencing formulas, however they are hardly used in practice, because of the 
advantage that central differencing has of being very accurate. Error bounding is not an issue for the 
diffusive fluxes. 
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6.4 EXAMPLES 
6.4.1 Stability region of a simple ODE 


For illustrating the above principles consider the simple ODE and its initial conditions as defined by 


equation (6.23) below: 
du(t) _ 0 ¢final 
a Au(t) fort e [r9, rf" ] (6.55) 


u(t®) = u? 


In the ODE defined by (6.54), A may be a complex number (i.e., à =A, + iA;). Introducing an equally 
spaced grid with (N + 1) nodes of discretization, such that: 


: ò A A pfinal - t? (6 56) 
peg +AR Ape. : 
i N 


and using an explicit discretization approach (Euler’s forward finite difference method), the solution of the 
ODE is given by: 


u"! = u” + At- À- u” = u”(1+ À- At), n=0,1,2,...,N (6.57) 
Equation (6.57) has a recursive character, hence: 

u” = u? (1+ A-At) = u? - An (6.58) 

where A is the amplification factor, a complex number of the form: 

A=1+A,-Att+i-A,-At (6.59) 
As previously shown, in order to obtain the region of stability the condition to be met is: 

IAI < 1) e IAP = (1+ ARY - (AA <1 (6.60) 
This condition represents the unit circle centred at Z = —1 in the complex plane Z. As long as À - At is 


positioned inside the circle the numerical method is stable. This shows that if A is a negative real number, 
stability condition requires that 


At € > (6.61) 
Such an expression shows that the method is conditionally stable. 


6.4.2 Convergence of an ODE: Emptying of a groundwater reservoir 

In order to study the accuracy of explicit and implicit schemes for ODEs the example of a cylindrical 
groundwater reservoir is considered. The volume balance of the reservoir (Figure 6.9) is given by a simple 
ODE, showing the variation in time of the reservoir water depth. 
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Figure 6.9 Cylindrical ground water reservoir. 


Assuming a constant reservoir surface area A, a constant porosity n and an outflow discharge Q defined 
as linear function of the water level h, gives the following ODE representation for the reservoir: 


dh 
2 -Q (6.62) 


Eliminating Q, which is a linear function of h, yields: 


dh o (6.63) 
dt 


where a is the linear reservoir coefficient. Equation (6.62) has the exact solution 


h = he (6.64) 


6.4.2.1 Explicit schemes 


An explicit scheme for equation (4.14), written on an equally spaced discretization (A — f) of the independent 
variable, gives: 


hj" = h" o At-h" = (1 — At) - h” (6.65) 


If the numerical solution (6.65) of the equation (6.63) is represented graphically for different (Ar), for 
a given initial condition of h = 15 m (see Figure 6.10), a comparison with the analytical solution (6.64) 
indicates that the solution decreases too slowly when At is large. Moreover this simple example shows 
that numerical methods are not always accurate, for example for a value of At of 22 days, failures in 
determining the accurate solution may be observed. This can be avoided by choosing the time step in a 
proper way. Analysis of the accuracy of the method is further done. 

Equation (6.65) is an accurate approximation of the equation (6.63) if it behaves in the same way as 
the function to which it is an approximation, and if it becomes a better approximation as the time-step is 
reduced. Starting from the first time-step of discretization, when f = 1^, equation (6.65) yields: 


h! = (1 — ot)? 


3 c: _ DL = E OP .— - 240 
h? = (1—aAtyh' = (0 «AD — aD? = (0— xA? h (6.66) 


h” = (1) At n? 
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Figure 6.10 Exact and explicit solutions for the reservoir equation. 


Relation of type (6.66) is analysed in behaviour as At — 0. Because At = t"/n, equation (6.66) can be 
written as: 


W= (2) h? (6.67) 


Using the binomial expansion relation (6.67) is: 


2 
h? = h? L + qt EIN + yr? - z| 8 € 1 (6.68) 


where A) is the number of ways of choosing r items from a set of n items. Order of choosing is not 


important. The number of possible combinations is: 


irr 6.69 
r) (n-r)r! (6.69) 
As n — ce, which is equivalent with At — 0, n becomes so big that n — 1 = n, n — 2 =n, and so on. As a 
result: 
n! n n” 
= n(n—1)(n-2)..—n and > — (6.70) 
(n — r)!r! r n! 
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Hence: 


2 v 
n 0 n n _ Ot pu _ t d p A i) n _% ES _ ot eri 
eoe ot] an sl JE "pl t+ T (671) 


at 


Taylor series for e^ 
Hence: 
hn + h9e-* (6.72) 


which is the exact solution of the ODE, that is, equation (6.64). 

The exact solution to the equation (6.63) decays with time. The value of h” at a given time step always 
smaller than the value at the previous time step. The decay is monotonic, i.e. the sign of h” does not 
change. Looking at equation (6.65), and taking into account the previous statement, if h”*! is less then 
h" yields: 


O<sl-aAts1 € OSATS (6.73) 


ZL 


Equation (6.73) shows that the numerical solution decays monotonically only for a limited range of 
values of At, as it could be noticed from Figure 6.10. 
Stability of numerical scheme (6.65) is analysed by applying the von Neumann principle: 


hy 
h” 


<1 © |-oAts1 (6.74) 


For this to be true, the two following conditions must be verified simultaneously: 
6.75 
Q ^t € 2 Rm 
-«^t € 0 


The first of the condition (6.74) is always verified because both the constant œ and the time step At are 
positive. The second condition means that the time step cannot be larger than a maximum value defined by 


2 
Atus = = (6.76) 


Hence there is a stability constraint attached to the time step. 

Consistency analysis aims to find out how accurately the discretised differential equation approximates 
the real differential equation. The analysis is done by using Taylor series expansions, expressing h”*! as a 
function of A": 

2 42 
dh Ať d'h (6.77) 


n+l — pn P raa 3 
hu =h TOI 2 PORA 
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The term O(AP) accounts for all the higher-order terms in the Taylor series expansion. Substituting the 
expression (6.77) into the numerical scheme (6.65), yields: 
dh At? d*h 


n paid! rp ae Ij = n (6.78) 
h” + At di + 2 d? + O(At^) = (1 — KAt)h 


Simplifying h” and dividing by At, yields: 


dh u At d'h 3 (6.79) 
“oo 8 gp t OA ) 


Comparing equation (6.79) with the original equation (6.63) shows that the discretised equation (6.79) 
is different from the original equation. The difference is the truncation error term, which becomes smaller 
when Aż decreases to zero. When Aż tends to zero, the term O(A?) decreases faster than the term Atd?h/dt 
hence the term in Af is the larger of the two. The truncation error is of first-order with respect to time 
step. The truncation error contains a second-order derivative d?h/d?, which is not present in the real 


equation, consequently, the discretised equation approximates the real equation up to the first order only. 
The discretization is first-order accurate with respect to h. 


6.4.2.2 Implicit schemes 


The implicit scheme for the reservoir equation (6.63) is: 


ya. P (6.80) 
1+ «^t 


Stability analysis of the implicit scheme gives 


h”! 
h” 


1 
1+ «^t 


<1 © |l+aAr21 (6.81) 


Condition expressed by (6.81) is always verified. Since œ > 0, condition (6.81) is valid for all At > 0 and 
At > 0 by definition. The scheme always decays and never oscillates around 0. It is unconditionally stable, 
hence there is no stability constraint for the implicit method. (Note: because a scheme is stable, does not 
make it accurate.) 

From physical point of view the reservoir becomes empty over time, hence h should decrease 
monotonically to 0, which is shown by the analytical solution. For the case of the explicit scheme Ah is 
invariably overestimated, to the extent that in some cases h is negative. In contrast, the implicit scheme, 
invariably underestimates the solution. 


6.4.3 PDE: Convergence analysis for Preissmann scheme applied 
to advection equation 


As seen there are a variety of numerical schemes that can be used to approximate the solution of the linear 
advection equation. One important scheme is the Preissmann scheme, also called the four point scheme or 
the box scheme. The four point scheme is often used for open channel flow modelling hence the study of its 
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convergence is an important aspect to be looked at. Recall from chapter 4, section 4.5.2.4, using notations 
from Figure 4.18 and equations 4.56 for the scheme that problem to be solved is: 


e + a =0 

ot Ox 

du ğ ut — ut tü-w) up*! — up (6.82a,b,c) 
Z .620,D,C 

ot E At At 

dul up -ut | jun! — unt! 

ax, Oe t 


The domain of computation of equation (6.82a) and the 4 point scheme is represented in Figures (4.19) 
and (4.18) of Chapter 4. The box scheme approximates the derivatives at a point P inside the discretization 
grid, using the values of the unknown function at grid points. There are in total M points of discretization. 
Advection velocity a is assumed constant and positive (i.e., a > 0), which implies that boundary condition 
is specified at x = 0, hence all values uj are known at all time levels n of computation. The two weighting 
parameters 0 and y ranges in the interval [0, 1]. 

Substitution of the derivatives discretization in advection equation yields: 

uji — un 


At 


n+l yn n y” n+l _ „n+l 
y a= yji —_" ZU peu y ADU |= 0i = nsn (6.83) 


At Ax 


In each of the (M—1) equations of system (6.83) there is one unknown, hence each equation of the system 
can be written in a general form as: 

ur = au, + Bur ty um (6.84) 
with values of the coefficients: 


y * (0 - 6)€r 


w (6.84a) 
l-y+0-Cr 
_1l-y -(1-0)Cr 6.84b 
B 1-yw48-Cr ia 
"m iiio 6.84c 
F7 1-yc8.Cr bd 
Courant number Cr is as previously defined: 
Cr- ^ (6.85) 


Ax 


Based on the linear system of equations (6.83) all values of the unknown function u(x, f) are computed 
at the time level (n + 1) provided that values at time level n are already known. Computations should be 
carried out in a consecutive order starting from the left boundary and advancing towards the right boundary. 
It is not possible to reverse the order of computation. Because of such constraints (e.g., computation order) 
the scheme is an implicit one. Consider the example from Section 6.2.2, a 10 km long canal with an initial 
water level of 2 m and the input stage hydrograph as defined by relations (6.26). Computation of the 
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solution is carried out over 10 hours. Results of solution for this problem, using numerical scheme (6.84) 


and different values for the weighting coefficients (0 and y), are shown in Figures 6.11, 6.12 and 6.13, at 
the discretization grid point x = 5 km. 


6- 


55. x=0 x=2500m 
` Boundary Condition 


Water Level (m) 


0 500 1000 1500 2000 2500 3000 3500 4000 4500 5000 
Time (sec) 


Figure 6.11 Solution of advection equation using the 4-point numerical scheme for 0 = 0.5 and y= 0.5. 
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Figure 6.12 Solution of advection equation using the 4-point numerical scheme for 0 2 0 and Cr- 0.5. 
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Figure 6.13 Solution of advection equation using the 4-point numerical scheme for y= 0 and Cr 0.5. 


Figures 6.11, 6.12, 6.13 shows that there are values of 0, y, Cr for which solutions are either oscillating 
or dampening. The physical phenomenon represented in the equation is just advection, hence the errors 
in solution are of numerical nature. Analysis of the error of the solution must give the constraints of 
application of the numerical scheme. 

As shown by relation (6.33) the error of the numerical scheme solution satisfy the equation of the 
applied numerical scheme, hence in the case of box scheme the error is: 


TET! = @-TE*, + B- ET? + y TET (6.86) 


(Note: Because Fourier series are used further for the analysis of the error, indices j are used for space 
discretization points, in order to avoid confusion with the imaginary unit i.) 
Using Fourier series to represent the solution error yields: 


A”! . ef PF = y. An . eh PUD + B . A” gii + y- A”! . ef PU-D (6.87) 
From equation (6.87) the amplification factor is computed as: 


And E p +Q- ei? 


A” ~ s 1 a y » ei? (6.88) 


The amplitude of the error should not increase, it is required that IAI « 1, hence: 


_ B+a-cosp—i-a-sing 


A : : (6.89) 
1-y-cosp+i-y-sing 
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Rearranging and using properties of complex numbers, one obtains 


p RU LT S (E CORNA A (6.90) 
1- 2y: cosp +y 1- 2y: cosp +y 
Stability condition yields: 
2 : 2 
Al = de Tae " (a+ B - y)sing «1 (6.91) 
1-2y-coso+y 1- 2y -coso + y? 


In inequality (6.91) pe [0, 7], hence the inequality is checked for the minimum and maximum values 


of @. 


If = z, inequality (6.91) is: 
e 2 
z) <1 (6.92) 
l-y 
equivalent with 
daf Sai (6.93) 
l-y 
Replacing the expression of coefficients œ, D, and y, yields: 
1< 2Cr 
— 2y — 2(1 — ~ 1-2y+20-C 
(a ey at j jd r (6.94) 
1- 2y +20 -Cr 2Cr E 
1-2y -20.Cr | 
If = 0, inequality (6.91) is: 
+a 
LOK 
= (6.95) 
equivalent with 
+a 
-1< < 
1< l-7 1 (6.96) 
Replacing the expression of coefficients œ, D, and y, yields: 
0 2+0-Cr 
1+2y ~ 1-2yw+6@-Cr (697) 
^1-2y *0.Cr ^ Ay + 0-Cr l 


1-2y +20- Cr ^ 
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Relations (6.94) and (6.97) are satisfied for: 


(6.98) 


2.20. 


1.60 + 


1.00 + 


Amplitude factor 


2.00 - 


Amplitude factor 


Figure 6.14 Amplification factor versus N for advection equation for Cr= 2: (a) different 0 values and 
y — 0.5; (b) different y values and 0 - 0.5. 
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Condition (6.98) ensures the denominator in relations (6.94) and (6.98) is not zero. Moreover the 
inequalities are satisfied for any value of Courant number, Cr, hence the scheme is unconditionally stable 
if (6.98) is valid. 

In the case that values of parameters 0 and y are zero the explicit FTBS scheme is obtained, which has 
different numerical properties than the box scheme, as shown earlier in the chapter. 

Further it is important to analyse the conditions for which box scheme does not introduce artificial 
dissipation or dispersion, because the analytical solution of the advection equation is a wave that has a 
specific amplitude and celerity which should not be changed by the numerical solution. The analysis is 
carried out, by examining the amplitude and phase errors. The amplitude errors are given by the modulus 
of the amplification factor, which for box scheme is: 


2 2 2 2 5 
B+ (a- B po) -a-y (a+ pysia 2) 
lAl = 5 + 3 (6.99) 
m 1 
1-2y- cos 27) Ly? 1-2y- cos 27) Ly? 


Recall that N is the ratio between the wave length and the size of the space discretization (Ax). 

Representing graphically the function A given by (6.99) for different values of Courant number (see 
Figure 6.14), it can be seen that for Cr = 1, and 0 = y= 0.5 the box scheme is free from dissipation, while 
for, 0 2 0.5 and y 2 0.5, the scheme increases the wave amplitude over time, which leads to instabilities. 
Small values of N and large Courant numbers creates damping, hence the number of discretization points 
over the wave length should be carefully selected. 


6.4.4 PDE: Convergence analysis for diffusion equation 


Analysis of convergence for diffusion equation using a FTCS scheme is analysed. The equations are: 


o? 


Ou u 

CAR oe 

Qu  u"-uw 

> EE a (6.100a, b, c) 
Ou uh, — 2u; tup, 

x P Ax 


where D is the diffusion coefficient. Substitution of the derivatives discretization in advection equation 
yields: 


DAt 
ut = u} +r- (ul, — 2u? + ut), withr = A (6.101) 


The error of the numerical scheme solution satisfy the equation of the applied numerical scheme, hence 
using Fourier series to represent the solution error yields: 


Atl. eS = AML + r (eh + ett — 2)) (6.102) 
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(Note: Because Fourier series are used further for the analysis of the error, indices j are used for space 
discretization points, in order to avoid confusion with the imaginary unit i.) 
From equation (6.102) the amplification factor is computed as: 


Vu 
"3 A-i-2r-(1-cosg) (6.103) 
The amplitude of the error should not increase, it is required that IAI « 1, hence: 
(6.104) 


r < 1/2 


which shows that the scheme is conditionally stable. 
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Chapter 7 


River system modelling and flood 
propagation 


7.1 INTRODUCTION 


Rivers are very important for mankind because they provide a source for drinking water and food 
production. Moreover they are important as a habitat for wildlife and for maintaining the equilibrium 
of the ecosystems formed along them. Apart from these natural benefits, rivers can also be used as 
means of transportation of goods and as recreation areas. Though sometimes rivers are causing floods, 
human settlements were always located in the river floodplains because of the social and economic 
benefits that the rivers are bringing (Knight, 2006; Novak, 1990). As societies developed their need for 
managing the river systems along which they live grew and changes started to be implemented in order 
to take advantage of the river's water availability or to protect the settlements against flood hazards. 
Communities have been developing and managing water resources in river basins in order to satisfy their 
demands for water quantity and quality or to cope with extreme events such as floods and droughts, by 
building dams to have access to water, irrigation channels to take water from rivers, dykes to protect the 
floodplains from flood risks, and so on. These developments were challenges to hydraulic engineers who 
needed to understand and evaluate the mechanism of flow in a river as well as the mechanism of flood 
formation. Stream flow predictions at particular locations in a river channel were done by engineers since 
a long time, and lately due to the advent of computers their tasks is facilitated by the availability of river 
system modelling tools, that are able to help in long and short term planning. The concept of River Basin 
Management (RBM) emerged and is accepted as a central element in all policy directives regarding 
rivers worldwide. The concept emerged as a holistic approach to river systems management and tries 
to eliminate the negative effect of the competitive interventions in a system due to the different water 
resources needs, such as unbalanced water resources exploitation between upstream and downstream of 
a river system. 

In the last four decades model codes have been developed for water resources and environmental 
analyses. These codes have been developed by academia, governmental and consultancy agencies and 
they vary in complexity from simple conceptual codes to complex ones based on physical processes. 
Majority of these codes are nowadays embedded in a tool, with complex graphical user interfaces. 
Further more, many of these tools are embedded in Decision Support Systems (DSSs), that are not only a 
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framework for providing data and information to decision makers, but they are also platforms for sharing 
and disseminating information to the public (Popescu et al., 2012). 

The DSS systems applied in RBM are using river system modelling, because it allows for understanding 
the behaviour of the river system as a response to some input variables such as discharges in an upstream part 
and/or conditions at the downstream boundary of a river. Moreover water quality analysis can be carried out 
if the hydrodynamic of flow in a river is known. A model of a river can help the decision making process by 
equipping the decision makers with a tool to objectively evaluate the effect of certain decisions. Mathematical 
(hydraulic) models for rivers, channel and pipe systems, emerged in 1980 and are the most widely used codes. 
They are used for predicting flows, travel times, and water level variations in rivers, channels and canal systems. 

Floodplains along the rivers are also important because they are the preferred location for living, 
however the effect of floods are loss of life, loss of property and adverse impacts on economic activities. 
Hence there is the need for accurate and reliable flood analysis and management of flood from both 
forecasting and mitigating points of view. Hydraulic models of river flow simulations became important 
instruments to address these issues. 

Development of river models started during 1970s. River models are in general part of river basin models, 
which were developed especially for planning purposes (Zagona et al., 2001; Jonoski & Popescu, 2004), 
both long and short-term. River basin models are used for water management purposes using different 
approaches and technologies. There are a number of models able to represent different hydrological 
processes, as well as river hydrodynamics. Examples of such systems that include sophisticated graphical 
user interfaces are MIKE SHE (DHI, 2012), HEC suite of models (USACE models) and SOBEK (Deltares, 
2010). There are many other examples of such systems. A comprehensive overview of existing river 
modelling systems is given by Horrit and Bates (2002) . 


7.2 RIVER SYSTEMS MODELLING 


Computational Hydraulics is a valuable tool for designer engineers because it allows for prediction of flow 
in a particular river and under particular boundary and initial conditions without making measurements 
in the field. The application of Computational Hydraulics to river related issues covers a wide range of 
problems such as flood alleviation, river rehabilitation (water quality and quantity), operation of reservoirs, 
and morphological predictions. 

The equations that govern the motion of flow can be solved analytically only for particular very simple 
problems, therefore the use of numerical methods for solving these equations is a good way forward. Rivers 
are having irregular topographies, which are leading to mixed regions of flows. This poses challenges to 
the numerical solutions of the flow equations, because they do need to solve these problems under a wide 
range of boundary and initial conditions. Moreover flow around and through structures (weirs, gates, etc.) 
needs to be addressed as well. All numerical approaches are subject to stability constraints, which restricts 
the selection of the time and space step and can extend the computational time. 

A number of numerical approaches exist to solve the flow equation processes. Corresponding to 
the approaches presented a number of mathematical models exists for solving flow in rivers. These 
are ranging from 1D Saint-Venant to 3D Navier-Stokes equations. All equations are based on the same 
physical principles and the selection of one model for a particular problem depends on the actual problem 
considered and the requirements on the solution. This chapter shows the numerical solution approach for 
the one dimensional Saint-Venant equations, which are simple as compared with the complex Navier- 
Stokes equations, however they are able to predict enough useful information for river engineers. 

Until recently the non-conservative form of the shallow water equations was used to compute flow patters, 
however this is not useful for problems that involve rapidly changing transitional flow, for which sharp 
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fronts and rapid changes in velocity occur. In this section of the book the conservative equations of flow are 
considered and solutions are shown. All demonstrations are made on rectangular grids because historically 
they have remained a popular choice and are the basis for a number of widely used algorithms (Sturm, 2001). 

For channels of irregular topographies, the flow can be described by the Saint-Venant equations, as 
shown in Chapter 2. The Saint-Venant equations for one-dimensional flow are: 


— continuity equation (CE): 


00 dh, 1 àQ 9h q a 


— momentum equation (ME): 


2 
90 4 2( aS} + ease AS, + 8A SS =0 (7.2) 
where Q is the discharge (m?/s); b the storage width at the water surface level (m); h water depth (m); A 
the cross-sectional area (m7); q the lateral discharge (m?/s); K the conveyance (m?/s); B the Boussinesq 
coefficient (2); g gravity (m/s?) and S, the bed slope of the channel (>. 

The system of equation formed by (7.1) and (7.2) has an analytical solution just for very few particular 
cases. These are usually solved using numerical approximation. There are many approaches to solve the 
system of equations however two of them, based on FDM, are the most used ones in practice, and are 
detailed further, the four point Preissmann scheme and the six-point Abbott Ionescu scheme. 

Prior to building any mathematical model and performing computations for a specific river, the geometry 
of the river needs to be known (Figure 7.1), by determining the course of it (the river bottom in a number 
of cross-sections. River systems are schematised into a simplified river network using a node-link structure 
(see Figure 7.1). The river network begins and ends with a node, and all nodes are interconnected. From the 
map representation of the river network, the extent of the model is selected by positioning the upstream and 
downstream boundaries. A computational space step Ax (usually constant) is selected as well. The obtained 
nodes represent the computational points on the grid and all relationship between the water depth A in these 
nodes and a number of water depth related parameters, such as width of water table (B), or cross-sectional 
area (A) are computed in all discretization points. Computational nodes do not usually correspond to the 
measured cross-sections and interpolation is used to determine values of parameters in computational nodes. 


7.2.1 Preissmann solution 


The Preissmann scheme for solving one-dimensional Saint-Venant equation was first introduced in 1961 
(Preissmann, 1961) in ‘Propagation of translatory waves in channels and rivers’. Preissmann scheme, 
also known as the four-point weighted implicit or box scheme is the most popular scheme used. Popular 
commercial software codes for solving river related problems, based on one dimensional Saint-Venant 
equations, such as HEC-RAS (Brunner, 2010) and DUFLOW (Clemmens et al., 1993) are using Preissmann 
scheme for solution of river flows. Preissmann scheme is suitable for relatively flat bed slopes. Dry-bed 
conditions are not allowed to occur during simulation in the codes using Preissmann scheme (Ogden & 
julien 2002). The scheme is applicable for both subcritical and supercritical flow, however the use of the 
scheme should not be extended to trans-critical flow regime (Chau, 1990; Meselhe & Holly, 1997). The 
limitation is due to the requirements on boundary conditions for different flow regimes. Abbott et al. 
(1991) shows that this issue can be solved by minimizing the effect of the convective acceleration terms 
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in momentum equations in order to ensure that flow remains in sub-critical regime. Hence the convective 
acceleration term is multiplied by a factor (Djordevic et al., 2004) which ensures that flow is always 
subcritical. 


River network 


e River link node 

Computational node 

o Cross section representative 
points 


Cross section points 


x 
m Boundary 


Cross section 


Datum 


Figure 7.1 River network, cross-sections and computational nodes. 


Discretization of the equations (7.1) and (7.2) is done on a computational domain (x, f), as it is represented 
in Figure 7.2. using the stencil from Figure 4.18 of Chapter 4 of the book. 


x,@ -unknown values at 
time level (n4 1) 
o -known value at 
time level n 
E -known value at 
time level (n4 1) 


Figure 7.2 Computational domain for the Saint-Venant equations for Preissmann scheme solution. 


The discretization grid considered is a structured one that computes values of Q and A in all points of 
the computational domain. Solution of the system of equations defined by (7.1) is solved for x € [0, L] and 
t 2 0. Applying the Preissmann scheme to the terms in equation (7.1) yields: 


oQ 4,90 dQ} 405 -Qr Qi - Q? 
d$ su dr 8 Ge $c (7.3) 


ox — 
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oh oh 


Z hii! — h; AAT 
ot ^ gi 


- y 25 i+] +(1 V) 1 i (1.4) 


oh 
*ü- 


i+] 


i 


Substituting (7.3) and (7.4) in equation (7.1) and grouping the unknown variables, Q?*',h?*', on the left 
hand side results in: 


Al,Q?*! + Ble! + CL Om! + DIA! = EI, (7.5) 


with values of the coefficients, in case q = 0, given by: 


-0 

Ala = <= (7.6a) 

b 

Bl = 07 V). (7.6b) 
0 

Clu = 4- (7.60) 
by 

"Li 7.6d 

Dl =F; (7.6d) 
1-0 b 

El, = (Of - Oja) + Geld — VOR + vtta (7.66) 


Applying the same approach (substituting (7.3) and (7.4) in equation (7.2)) to momentum equation, 
yields: 


— Discretization of first term of momentum equation: 


9Q Qi — Ohi, Qr" — QF 73 
w V At td WG et 


— Discretization of second term of momentum equation: 


TORORO] 


Considering the following linearization for the non linear term Q? (Abbot & Basco, 1989): 


Q? = Q" x Qn 
and 
n n+l 
A. = Anti = Aa + Am 
i+] ^ ^i = 2 
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the second term of the momentum equation is: 


a (0?) _ g[guxom _ oF x or 
Ox A Ax (Ani) Axl An?) 


— Discretization of third term of momentum equation: 


oh oh 
gu om = gAn? Ç Pm 


LN 


n+l 


hu —-qgmnm hp — hp 
IL AS *0-674 | 


— Discretization of fourth term of momentum equation: 


gASo = BANi” So 
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(721b) 


(70) 


(77d) 


— Discretization of fifth term of momentum equation is done using the coefficient y because Q and K 


are varying over time: 


PCR AR 

i+] i 
hence: 

n n+l à 1 n ntl 
gA aS- = g WA, 5 ; Qa Qu + 1 )A; 2 Q; = 


al ( 7 
ie] 


Grouping the unknown variables, Q"*! , "*!, on the left hand side results in: 


42,0"! + B2,h! + C2,Q!:! + D2, het! = E2, 


with coefficients: 


a l-y sa +a yan E 
1 At Ax(A? ) (K+?) 
[ 0 

B2, = | -8A 2l 
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(79a) 


(7.9b) 
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y Qi n+1/2 Qi 
Cha =| + BL + gyAT n (199) 
WU DAE T Axe = (Kr? ) 
1 
een ey i. (7.9d) 
i+] gy i+] Ax h 
i 1-90 
EA, = OF +E On + gAn EIC - hf) + gARIPSo (799 


The pair of equations (7.5) and (7.8) are written for each point of the discretization grid, that is for 
i=1,2...,M-—1, to obtain a system of 2(M — 1) equations that has 2M unknowns. The system is completed 
with two more relations at the boundary of the domain: 


— for node i= 1: 
bc(0) - h+! + (1 — bc(0)) - Q'*! = bc(0) - h (£1) + (1 — be(0)) - Qu (P) (7.10a) 
— for node i = M: 


bc(L) - hy! + Q- be(L)) - Qs" 
= bc(L) - hy (t"*!) + (1 — be(D)) - Qy C) (7.10b) 


In equations (7.10) bc(x) is a coefficient which takes the value O or 1 depending on the type of the 
specified boundary condition, that is coefficient value 1 means water depth is specified as boundary 
condition and value zero means a discharge boundary condition. 

A number of 2M equations are obtained for a number of 2M unknowns. Such a system can be solved 
using matrix calculations or in an iterative method by using the so-called Double Sweep Algorithm. 

Solution of the linear system of equation is given by the solution of the bandwidth matrix: 


| Left BC Q, | | Valuefrom BC | 
Al, Bl Cl, Di, xl h El, 
A2, B2 C2, D2, tae Q, E2, 
Al, Bl, CL Dl, -- h, El, 
A2, B2, C2, D2, + E2, 
© Aly Bly Clu Diy Ely. 
| A2y B2,4 C244. D2y4 Qu E2y-1 
L Right BC || hy | | Valuefrom BC | 
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An alternative solution to matric calculations is the one given by the Double Sweep Algorithm, which 
takes into account that at each point of discretization grid a relationship can be written between the 
unknowns, Q"*! and //*!, as well as between , Q”*!,A?*!, in the form: 


Qn 
i 


= Eh" +G, (7.11a) 


hU = HN + hA + J, SSH) 


Coefficients H, I, J, F, G are derived by replacing relations (7.11) into the system of equations formed 
by equations (7.5) and (7.8): 


zl 
«ifi: 
El; A6. 
j= a (7.120) 
= DAX(ALF, + Bl) - D1,(A2,F, + B2;) _ -1(A2;F, + B2;) - D2, im 
^ 7 CL(G;F, + B) - CAL, + Bl) | H(A2F, + B2) + C2, i 
E2; — A2, F,J; — A2;G; — B2;J; 
Gin, = i baa DS ivi (7126) 


H,(A2, F, + B2,) + C2, 


The step by step procedure to solve the system of equations for time level (n + 1) using the double 
sweep algorithm is illustrated in Figure 7.3 and follows the steps: 


(I) 
Q) 


(3) 


(4) 
6) 
(6) 


(7) 


Initial data in terms of water level and discharge are specified for all grid points from point 1 to 
point M of discretization grid, at time level n. 

In order to start the computation of Q and h at time level (n + 1) the assumption that all values 
of water level and discharge at time level (n + 1) are the same as the values of water level and 
discharge at time level n is made. 

Computation of all A1,B1,C1,D1,E1, A2,B2,C2, D2, E2; coefficients at level (n + 1) is done 
using formulas (7.8a—e). Coefficients are computed based on available data at time level n and 
(n+ 1). 

Boundary data at the upstream end is provided in terms of discharge or water level in grid 
discretization point 1. Coefficient bc(0) (see equation (7.11a)) is set to the appropriate choice. 
Values of Q and h at time level (n + 1), at the left boundary (upstream) are used to compute 
coefficients F and G in point 1 of the discretization grid (i.e., F, and Cj). 

Forward sweep: Using recurrence relations (7.12d—e), values of coefficients F, and G; are 
computed by sweeping across the grid (from point | to point M). Coefficients H,, I, and J; are 
calculated using relations (7.12a—c). All values are computed at discretization points of the grid, 
as opposed to coefficients A,B,C, D, E, which are computed in grid cells (see Figure 7.3). 

Based on downstream boundary condition, which may be either discharge or water level the other 
unknown (water level or discharge, respectively) is computed using relations (7.11a). 
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(8 Backward sweep: Using previously computed coefficients F, G, H, I; and J, discharge and water 
depth at time level (n 4 1) are calculated by sweeping back across the grid (from point M to point 
1) calculating discharge and water depth at time level (n + 1) using relations (7.12a—b). Discharge 
and water depth are determined at time level (n + 1) for all points. 
(9) Computed values for water level and discharge, obtained in step 8, are compared with the values 
from previous iteration. If they are not in accordance with a comparison criterion step 10 is applied. 
(10) Iterate: Computation is repeated starting at step 3, using the values obtained in step 8 to re-compute 
all the coefficients in a forward sweep. 
(11) Iterations are repeated until a certain condition is fulfilled (difference between two computed 
values, maximum number of iterations, etc.). 


A Compute h,Q 
t » 
n+l n+l 
Cohpute AB, C, »| E, F, cn 
Grid 
Uo i cell i 
n n, ul 
] il d il M ; up 


Figure 7.3 Double sweep algorithm. 


7.2.2 Abbott-lonescu solution 


A second FDM approach for modelling of unsteady river flow is the Abbot-Ionescu scheme, which is 
implemented by Danish Hydraulic Institute in their MIKE 11 river flow model component of Mike Zero. 

The scheme is based, as shown in Chapter 2, on a staggered grid (alternating grid-points). Based on 
this approach the solution is sought as Q and h values at alternate points of the discretization grid (see 
Figure 7.4). 


x,@ - unknown values at 


time level (n+1) 


n+l |Q h OQ h |Q |h 


O -known value at 
time level n 


Q h 


h |Q |h gi -known value at 


time level (n+ 1) 


X 
Figure 7.4 Computational domain for the Saint-Venant equations for Abbott-lonescu scheme solution. 
Figure 7.4. shows the one dimensional (x, f) computational domain for which the method is applied 


to solve the Saint-Venant system of equations. The Abbot-Ionescu scheme has been developed to solve 
subcritical river flow, hence one boundary condition has to be defined upstream and one downstream. For 
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showing the principle of the method the boundary condition upstream has been selected to be inflow (the 
most common condition) and the downstream one has been selected to be water depth, which means that 
the time series of inflow Q(x — 0, f) and A(M - Ax, f) are known before the computation of the solution, 
in the entire domain, starts to be determined. The number of discretization nodes along the x axis is M 
starting with 1 as the first point. An odd number of points are required for the discretization. 

The staggered grid definition and the boundary conditions imposes that at all time levels, in all 
even points of the discretization grid (1 = 0, 2, 4... , M — I) a Q value is computed; and in all odd ones 
(i= 1,3, 5,..., M), ah value is computed. 

In order to demonstrate the method it is considered that all values of discharge and water depth, at time 
level n, Q” and h”, respectively, are known; and they are the basis for computing the values at time level 
(n+ 1) where Q"*! and /*! are not yet known. 

The Abbot-Ionescu formulas, as described in Chapter 4, equations (4.61a—c), applied to continuity 
equation (7.1), by considering a weighting factor 0 — 0.5, yields: 


oQ 1 (Qr Q0 un -Q9u 
ara Bele 8 2.Ax (7.13a) 
Oh Ql (h — pe (be + bP) (noe 
VES m. 2 ow] i|- E pic i 7.13b 
j ot b; | At | 2 At ( ) 
m, qt tq dia + dia 
157A | z AER OA (7.130) 


Inserting equations (7.14a—c), in (7.1) and grouping all unknown terms, Q”*' and A""!, at time level 
(n + 1) on the left hand side of the equation, it yields: 


Al; «Qi * Bl; z. * Cl; «Qui = Di; (7.14) 


where the coefficients of Q"*' and A"?! are: 


At 
e 7.15 
dice x (7.15a) 
1 
Bl, = 5 (Br + BP) (7.15b) 
Cl, = —Al, (7.150) 
DI, = Bl; - h? + Al; - (Qf - QA.) + 4 (af att + aha + a) (laa) 


Equation (7.14) is linear in of Q”""! and h”*! and is valid for i = 1, 3, 5,..., M — 2. The equation is space- 
centred around a computational point where A is the unknown value. Because the weighting factor is 0.5 
the equation is also centered in time, between the time levels n and (n + 1). 
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Applying the Abbot-Ionescu formula of discretization for the derivatives of momentum equation, 


yields: 
oh D Qr Q” (7.16a) 
Ox At 


Əh 1 |O- Or Or - on, 


Ox 2- Ax nm "n (7.16D) 
il i-l 
dh my 1 (meth a — 
Bud SUE Y DES + 2. Ax (7.16c) 
up 
g-A-S, =g A, 2S, (7.16d) 
2- |o] Q -lQ "5 Or" - O; 
pa qeu ode A Cn (1169) 
M? - R3 - A? (Ky? 


where K = SM - R?? is the conveyance of the cross-section and SM is the Strickler coefficient. 

Equations (7.16a—e) are linear, because the term in Q, which was at power two is expressed as a product 
between a known value at time level n and an unknown value at time level (n 4 1). 

Similarly, using relations (7.16), for momentum equation, yields: 


A2, - htt! + B2, Q?" + C2, - het! = D2, (7.17) 


Coefficients of equation (7.17) are: 


1 
A2; =-g: A P iu (7.18a) 
n+ Q? 
B2, 210 € AL ? n (7.185) 
K*) ? 
C2, = -A2, (7.180) 
At nil, n nil , gn 
D2, = Q” , Qu Qa Qi Qr (7.18d) 
2. a Ax nts nt 
Ais Ai, 
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Equations (7.14) and (7.17) form together a linear system of equations of the form: 


[ Left BC Q, |_| Value from BC | 
Al. Al, C1 = h, DI, 
A2, B2, C2 pe Q, D2, 
Al, Bl, Cl, is h DL, 
A2, B2, C2, ee D2, 
- Alua Blua Clua . Dlui 
A2y- B2y.~  C2y1 Ou D2y-1 
L Right BC || Ay, | | Value from BC | 


The above system of equations can be solved using the classical approach of matrix algebra from 
calculus or the most used method, the Double Sweep Algorithm, as it was demonstrated previously 
for Preissmann scheme which uses the physical properties between Q and h to simplify the way the 
computation is done. As presented in the previous section yhe method of solution is based on continuous 
replacement from one equation to the next. In equation 1 from the linear system (7.14) the term Q, 
can be transferred in the right hand side in order to obtain an expression for h, as a function of Q, and 
other known terms. Further the expression of h, is replaced in the second equation and then into the 
following and so on. Eventually the last equation of the linear system contains one unknown value, the 
Qy. .. The value of the unknown is calculated and a backward insertion trough the equations can be 
started in order to determine all the unknowns. Computation is carried on several times, until a stop 
criterion is fulfilled. Moreover it needs to be noted that the coefficients are expressed using values of 
Q and h in discretization points where these are not computed, because of the staggered nature of the 
grid. This issue is overcome by using interpolated values for flow in grid-points where a water depth 
is computed and vice versa by using interpolated water depth in Q grid-points. The interpolation can 
be linear or parabolic, and can be calculated as an intermediate value or as a final step in the iteration 
procedure. 

As shown in Chapter 6 the solution accuracy is influenced by the choice of Ax and Ar. Given nowadays 
computer power the choice of Ax (implicitly of number of discretization points M) is not a problem 
anymore. Based on several application results as well as theoretical analysis conducted by Wright and 
Crossato (2013) concluded that the Abbott-Ionescu-Scheme gives numerical stable results without 
numerical diffusion. Chintu (1986) mention that the Abbott-Ionescu scheme works properly with Courant 
number values close to 10 due to the implicit nature of the scheme. Compared with an explicit numerical 
scheme, where the values of Q and A in a point P of the computational domain depends on the values of 
3 points from the previous time level, within an implicit scheme as is the Abbott-Ionescu scheme these 
values depend on 20 points, located to the left and the right of the P point, in the previous time level (see 
Figure 7.5). The Courant number can be considered as a measure of the extent of the flow regime from 
time level n to time level (n + 1). The requirement of Cr « 10 for the Abbott-Ionescu scheme is equivalent 


Downloaded from https://iwaponline.com/ebooks/book-pdf/650788/wio9781780400457.pdf 
bv IWA Publishina. publications@iwap.co.uk 


River system modelling and flood propagation 137 


to determine the values for Ax and Ar as: 


" > "T (7.19) 
Ah 
Accu g-h| < 10 (7.20) 
P 
n+l 
n 
| | | 
US 
" 2Ax i 
20Ax 


Figure 7.5 Domain of dependence for a point P in the computational domain, when Abbott-lonescu 
scheme is used for solution of Saint-Venant equations. 


The values of Ax and Aż as defined by equations (7.19—7.20) have to be correlated with the variations 
in initial and/or boundary conditions and of the tributaries of a modelled river. The choice of Ax and At 
should make a good description of all variations in times series as well as longitudinal profiles. 


7.2.3 Initial and boundary conditions 


Modelling river channels give good results if appropriate initial and boundary conditions are set to the 
model. Initial conditions are easy to set, and they usually come with the knowledge about the modelled 
system. These are given in terms of water levels or discharge values. 

Boundary conditions that are usually set at the upstream end of a river are water level or inflow variations 
in time. A rating curve should never be specified as a boundary condition at the upstream end of a model. 
Boundary conditions at the downstream end of a model are water levels, outflow hydrographs and rating 
curves. Cunge et al. (1980) advice that the same type of boundary conditions (Q or /) at both river ends is not 
recommended, because solution strongly depends of the initial conditions. The most used boundary conditions 
at the downstream end should be measured Q and h as variations in time. If discharge and/or water levels 
variations in time are not available at a downstream end of a river reach a rating curve (Q — h relationship) is 
used. However, such an approach may give problems, because it is represented as a single valued function. 
In practice a rating curve has hysteresis for unsteady flow conditions. In order to eliminate the influence of a 
rating curve Singh ef al. (1996) suggests that models should be artificially extended in the downstream. 


7.2.4 River networks 


It is important to mention that in case lateral inflow is considered in Saint-Venant equations, then relation 
(7.14c) takes into account the possible variations in space and time by considering the values of q from 
space point i to space point (i + 1) and values at time levels n and (n + 1). Notice that lateral inflow variation 
in time and space has to be known when applying Preissmann or Abbott-Ionescu schemes. If the model has 
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a network of river branches that are joining together, than a tributary is a lateral flow to the main branch. 
Figure 7.6 gives the possible variation of q in space at each time step. The flow values coming from a 
tributary at time level n (f" = n - Af) are added to the lateral inflow term. In case that Abbot-Ionescu scheme 
is used for the solution of the Saint-Venant equations due to the staggered nature of the discretization, grid 
tributaries have to be associated with odd i-points and they should be dispersed over a distance 2 - Ax. This 
approach requires careful selection of Ax. 


(a) Longitudibal profile and lateral inflow hydrographs 
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(b) Longitudinal profile and lateral flow at time level n 
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Figure 7.6 (a) Lateral inflow variation for network river systems and (b) tributaries handling in the Abbot- 
lonescu solution approach. 
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In case a river is formed from a network of channels the same procedure of solution is applied as for 
a single river reach, provided that boundary conditions are specified at all ends of the pendant branches 
that is in points 1, 71 and 101, as per Figure 7.1. A system of algebraic equations that need to be solved is 
obtained after adding additional equations for the channel junctions (see Figure 7.6): 


Q = Q; +Q; 
h, =h, (1.21) 
h; = h 


First equation in conditions (7.22) represents continuity equation, whereas the later two are simplified 
forms of energy equation written at the junction, under the assumption that local losses are neglected. 
Forsius and Huttula (1982) show that if flow velocities are too high, such that velocity head becomes 
significant, the equal water level condition in (7.22) should be replaced with an equal energy level 
compatibility equation. 

The corresponding matrix of the solution system for the Saint-Venant equations is represented in Figure 
7.8, for the river network of Figure 7.7. 


— 9 Direction of flow 


4Ax 


Channel C 


Figure 7.7 Channel junctions. 


The double sweep method can also be applied to network river systems as long as certain computational 
order is respected (Wu, 2008). For the channel network of Figure 7.7, the computational order is as 
follows; the forward sweep computes the recurrence coefficients for channel A starting from 1, where 
boundary condition is known, and is carried on down to junction point j. A similar sweep is carried on 
channel B from first to last point. At the junction the three cross sections should be located very close 
to each other and according to relation (7.21) it is assumed that the water levels at the cross sections 
are equal and discharge at the downstream cross section is equal to the sum of discharges in upstream 
cross sections of channel A and channel B, hence initial values for coefficients F and G for channel 3 
are computed. The forward sweep is further carried out from point i to point 3 of channel C and the 
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recurrence coefficients are computed. Backward sweep start from the most downstream point of channel 
C (point 3), where boundary condition are known. Values for discharge and water level are computed 
first along channel C back to the junction. From the junction the backward sweep algorithm is carried on 
further along channel A and B. 


Figure 7.8 Structure of a matrix for the river network system represented in Figure 7.7. 


7.3 MODELLING FLOODS 


Flood are phenomena characterised by high discharges and/or water levels, which may cause inundation 
of floodplains (the land adjacent to rivers) or to terrain close to water bodies. In general flooding happen 
when rivers are unable to route the amount of water coming to the river due to runoff. Because a river has 
no conveyance capacity leads to overtopping of river banks. Floods causes disruption of activities and loss 
of lives, hence they are given a lot of attention in water engineering practice. Flood disasters account for 
about a third of all natural disasters (by number and economic losses). 

Flood events may be caused by intense or long-lasting rainfall, snowmelt, dam or dyke break, ice jams, 
high tides, storm surges or by operation of reservoirs. 

One dimensional models for Saint-Venant applications have been widely applied in simulating flood 
routing and unsteady flow in river networks (Hicks et al., 2005; Gichamo et al., 2012; Moya Gomez et al., 
2013). Though there are advances in the two (2D) and three dimensional (3D) hydrodynamic modeling of 
flow in rivers, models based on one dimensional Saint-Venant equations still are the most popular choice 
for solving large-scale river engineering problems. The choice is due to their reduced data requirements 
compared to those for 2D and 3D models, as well as reduced computational cost. 

Special insights in flood propagation in river channels are gained by using Price’s (1985) analysis of 
the Saint-Venant equations (7.1), which can be rearranged in the form: 


oQ Q|1 dK 1 db|( ðQ K? yg E 
ar b E dh 2b A ox «| 2bQ Ox) — a (7.23) 
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Equation (7.23) is obtained by neglecting the inertial terms (the local acceleration and convection 
of momentum), and has the form of an advection-diffusion equation with advection and diffusion 
terms: 


_Q[idK  ! db 
«o, E dK 2b z| NE 
and 
K? 
DO.h) = 35 (724b) 


Coefficients given by relations (7.24) show that when a flow disturbance occurs on a river channel, 
a process of translation/advection takes place along with a corresponding diffusion. The phenomenon 
is known as attenuation of peak flow. The speed of the wave traveling downstream is given by the 
advection speed (7.24a). Note that equation (7.23) has Q as the primary dependent variable within the 
differentials, however because Q is dependent on A (or A), a continuity equation is required in order to 
solve it. 

An more general representation of equation (7.23) given by Price (1985) is: 


9Q 9 E 2) = Og — £C, = Qd J+ Ole?) (7.25) 


+ + 
g ^u ey c3 Ox 2qAs, 


where 
h 
gud. (7.26a) 
So 
_ @ 
y= C (7.26b) 
Q |, yB(. QY 
ay = 2s B I "m (a 2) | (7.260) 


Equation (7.25) is the basic flood routing equation, valid for any Froude number. It is also known 
as the Muskingum-Cunge equation. Muskingum is the name of a river in North America and Cunge is 
the scientist who first established the relationship between the Muskingum method and the Saint-Venant 
equations (Cunge, 1969). Based on equation (7.25) attenuation parameter can be determined for any river 
for specific discharge values. An example of a typical such curve is given in Figure 7.9. 

Flood attenuation studies (Anderson et al., 2006; Acreman et al., 2003; Castro Gama et al., 2014) 
concluded that the peak discharge attenuation is affected by floodplain (width and roughness) as well as 
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by bed slope. Hence floodplain storage is an important process to take into consideration while predicting 
the peak flood attenuation. One-dimensional models should be capable to account for floodplain storage 
during flood events, hence, considering a coupled 1D/2D model in flood simulations is important. An 
extensive overview of flood inundation models is given by Horrit and Bates (2003). 


aQ) 


M 


Bank full Q 


Figure 7.9 Typical flood attenuation curves. 


7.4 RIVER ROUTING EXAMPLE 


In order to show the capability of a river to route flow the example of a simple rectangular canal is 
chosen. A flood hydrograph is routed downstream on the 200 km long channel in different conditions. 
Different widths of the rectangular cross-section are considered, as well as different slopes and 
roughness. HEC-RAS is used (i.e., Preissmann scheme) to compute the different cases. The hydrograph 
that is routed downstream the considered canal is shown in Figure 7.10, along with the schematic of the 
canal geometry. 


Q 


Boundary 
condition 
upstream 


L= 200 km b 


Figure 7.10 Canal geometry and boundary conditions. 


Results of routing are shown in Figures 7.11, 7.12, 7.13. Figure 7.11 shows the effect of roughness in 
flood routing. Figure 7.12 shows the effect of the slope on river floods. Steep slope route quicker the flood 
peak, however the attenuation is smaller. Figure 7.13 shows the effect of cross-section width on river 
routing. The wider the river cross-section the higher the attenuation of the flood peak, however longer the 
time for the flood to be routed out of the river system. 
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Figure 7.11 Effect of roughness on routing a flood wave in a river; (a) n 2 0.025; (b) n 2 0.030. 
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Figure 7.12 Effect of slope on routing a flood wave in a river; (a) S = 0.0003; (b) S = 0.00001. 


Downloaded from https://iwaponline.com/ebooks/book-pdf/650788/wio9781780400457.pdf 
bv IWA Publishina. publications@iwap.co.uk 


River system modelling and flood propagation 145 


(a) Upstream & Downstream hydrographs U/S & D/S rating— curves 
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Figure 7.13 Effect of cross-section width on routing a flood wave in a river; (a) 50 m width; (b) 100 m width. 
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Chapter 8 
Water quality modelling 


8.1 INTRODUCTION 


Fresh water is very important for development because it is a source of energy, a mean for transportation 
and a habitat for aquatic life. Availability of good quality water (fresh water) is very important for economic 
activities and growth, however as population and economies grow, so does the amount of pollutants that 
are discarded in water bodies (Orlob, 1983). When levels for water quality are not met, treatment must be 
carried out, which brings extra costs to water users. 

The issue of water quality is not something new, as Heathcote (2009) notes Romans issued the following 
Roman decree: 


‘Ne quis aquam oletato dolo malo ubi publice saliet si quis oletarit sestertiorum X mila multa esto, 


which reads 


‘It is forbidden to pollute the public water supply: any deliberate offender shall be punished by a fine of 10,000 
sesterces'. Year 1989 was the 2000th anniversary of the issuing of this decree. 


Countries around the world defined standards for water quality and make plans for reaching good 
ecological status of their water resources. A good example of such efforts is the European Water Framework 
Directive, which is a major policy in the area of water resources management in European Union (EU). All 
countries of EU follows the policy, which sets clear requirements regarding attitudes and approaches to 
water resources at both national and regional levels. The directive's main goal is sustainable use of water 
resources with a dominant ‘environmental objective’ of preventing the deterioration of the status of all 
EU waters. 

In this context of continuous need to maintain and improve water quality scientists and managers 
are continuously trying to come with solutions for providing understanding of freshwater ecological 
systems and predictions of the water quality responses to human interventions, as well as evolution of 
these responses in time. Hence water quality models are important tools assisting in the analysis regarding 
pollution of water and ultimately in the identification of the state of the environment. Such tools need to be 
able to show the impact of humans activities on the status of a water body (i.e., lake, river, etc.). 
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Water quality modelling aims to simulate behaviour of variables related to quality of a water body. 
Variables include characteristics of the water body, such as temperature, solar radiation, wind speed, as 
well as pollutants existent in the modelled water body. Models of water quality describe physical processes 
in terms of hydrodynamics; and are complemented by a description of chemical and biological processes 
that control the transport and transformation of the variables. Pollution elements are defined as point- 
source and nonpoint source loads. Developed water quality models range from very general, applicable to 
all kind of water bodies; to very specific ones, that are applicable to a particular problem type. Hence each 
model may have its own characteristics and requirement for data input. Water quality models are known 
not to be able to accurately predict what happens, however they are still valuable to understand processes 
and estimate relative changes due to management practices and policies. Hence even the simpler water 
quality model is a useful tool to be used. 

First water quality model, the so-called S-P model, was defined by Streeter and Phelps (1925) for Ohio 
river, in United States. The S-P model is steady state one-dimensional determining the oxygen balance and 
the decay of bio-chemical oxygen demand (BOD). Since 1925 models have evolved from steady state models 
to dynamic ones and from point source pollution to non-point pollution sources (Jung et al., 2011), as it is the 
example of EFDC model developed by Virginia Institute of Marine Science, which is suitable to be used for 
water quality simulation in rivers, lakes, reservoirs, estuaries, and wetlands, for one-, two-, or three-dimensions 
(EPA, 1999). Table 8.1 presents the main available models and their characteristics, as they evolved in time. 


Table 8.1 Water Quality Models. 


Model Year of development Main characteristics 


Streeter-Phelps 1925 — One-dimensional 
— Steady state 
— Oxigen balance 
QUAL 1970 — One-dimensional or dynamic 
— Non-point source pollution 
— Suitable for dendritic rivers 
WASP 1983 — One-, two-, or three-dimensional models 
DHI models 1990 — One-, two-, or three-dimensional models; 
— Non-point source pollution; 
— Suitable for rivers, lakes, estuaries, and reservoirs; 
BASINS 1993 — One-, two-, or three-dimensional models; 
— Non-point source pollution; 
— Suitable for catchments; 
EFDC 1997 — One. two-, or three-dimensional models; 
— Non-point source pollution; 
— Suitable for rivers, lakes, estuaries, and reservoirs 


An important aspect of water quality modelling is the availability of monitoring data, which decide 
the appropriate model to be used for the complexity of the situation (Louks & Beek, 2005). Depending on 
the availability of monitored water quality data it is advisable to begin with simple models and gradually 
extend to additional complexity as data becomes available due to data collection campaigns. Prediction 
of hydrodynamics of the water body has long been the recognised as an important component of any 
water quality study (Martin & McCutcheon, 1999). It is well known that river have good collection data 
for hydrodynamic parameters, however most water quality monitoring programmes for water bodies with 
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stagnant waters, detailed hydrodynamic measurement is hardly included neither is prediction of water 
motion, the emphasis is mainly on general physical-chemical parameters due to their ease of measurement. 
In case of lakes water quality parameters collected during spring and fall overturn represent the best 
and most uniform water quality conditions and are most valuable for calibrating and validating models. 
However making model predictions based on one or two samples is not enough. More extensive sampling 
is suggested because it provides additional information. 

Present chapter presents the main equations that are used to define water quality models. First the main 
processes described by a water quality model are presented, followed by specific particularities of processes 
for rivers and lakes. Models used for lakes and reservoirs are substantially different from the ones used for 
rivers. Usually water movements are very slow in lakes and reservoirs, hence the hydrodynamics is driven 
by wind; heat exchange due to solar radiation; and inflows and outflows. Phenomena that are normally 
absent or negligible in river flow are very important in non-flowing waters (‘standing’ waters) such as 
lakes and wetlands. For example in lakes the phenomena of temperature stratification appears, which can 
influence the vertical turbulent mixing; or different density stratification of water causes advection. These 
phenomena are known to create a particular thermo-hydrodynamics state in a lake. In rivers the main 
driving force is gravity and bed friction (as resisting force), which generates perceptible velocities and 
turbulent mixing. Moreover sediment transport and deposition is different in rivers and lakes, therefore 
their effect on water quality is different. Usually in lakes and reservoirs a phenomenon of slow deposition 
of fine sediments takes place and the bottom of the lake can release substances, as effect of bacterial 
decompositions and chemical processes that are taking place. In rivers such processes are negligible as 
compared with the advective transport of substances in the water. 


8.2 PROCESSES DESCRIBED IN WATER QUALITY MODELS 


Water quality models are applied to different types of water bodies; rivers, estuaries, lakes, reservoirs, 
coastal waters and seas or oceans. Models are usually used to describe the main water quality processes, 
and require input in form of hydrological time series of inflow and outflows, as well as pollutant loading 
(pollutant concentrations). 

Water quality models are a set of equations that describe the advective and dispersive transport of 
pollutant constituents within the water body, based on the hydrodynamic conditions; as well as a set of 
equations describing physical, chemical and biological reactions among constituents. 

Advective transport is dominant in river, where flow takes place, while dispersion is predominant 
transport phenomenon in estuaries, which are under tidal influence. Lakes and reservoirs water quality is 
more difficult to predict than the one in rivers because of the wind direction, which is random in nature. 
Currents formed due to wind influence affect surface mixing and generate stratification of the water. 

The basic principle of water quality models is application of mass balance. Space computational domain, 
is formed by one-, two-, or three-dimensional system of the water body (i.e., river, lake, estuary, a zone 
in the sea or the ocean). Space and time domain is divided, using a discretization grid, in computational 
cells (volume elements) in space and in discrete time intervals Ar. Water quality is computed over time 
for a computational cell volume by writing the mass balance for each water quality constituent. Flow 
characteristics in the control volume, over the time interval Af, is pre-computed as a hydrodynamic solution 
carried out for the given inflow and outflow conditions. 

For each computational cell and each time step At, the mass balance of a substance can be defined as a 
sequence of step by step computations of: 


* changes due to transport in and out of the computational cell (TR), which can have two components; 
and advective term (flowing water) and a dispersive term (pollutant concentration differences); 
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* changes due to physical or (bio)chemical processes (P) occurring within the computational cell; and 
* changes attributed to sources/discharges to or from the computational cell (S). 


Mass balance over a time step interval Ar, for computational cell i at time level t + At can be expressed 
as: 


t+At t+At t+At 


M(t + At) = M(t) + J Mrpdt + J Mpdt + Í M;dt (8.1) 


In equation (8.1) M is the mass in computational cell i, M, is the mass change in computational cell i, 
due to component F, which is transport (TR); (bio)chemical or biological processes (P); and sources (S). 

Equation (8.1) of the mass balance for computational cell i at time level (n + 1) can be further expressed 
as (Louks & van Beek, 1983): 


MP" = M? +At- AM: TA AM, + At: AM, (8.2) 
At Jor At Jp At Js 


where M” is the mass in computational cell i at time level n and M! is correspondingly the mass in 
computational cell i at time level (n + 1). The term | AM, | represent the changes in computational cell i, 
due to component F. A 

Transport due to dispersion is important in lakes where water column is stratified, and it is less visible in 
rivers. Dispersion is not the same as diffusion. Figure 8.1 shows schematically possible transport elements 
in lakes and rivers. 


(a) 
water air exchan ges 
"ANsport/mixin (hori 
e | 8(orizontal & verti 1) 
ca ) 
sedimentatio 
n i 
sediment res i 
uspe: 
sediment — m 
"Transport processes in a river 
(b) 


water air interface 


surface water 


inflow horizontal transport/mixing 
€ » outflow 
g g a 
" m. E 5 
ground water E ^ z S., 
. m. en E B 
inflow E N- | PE 
i iz 9: 
9 MES EM 


seepage 
sediment water interface 


sedimentation 


Transport processes in a lake 


Figure 8.1 Schematization of transport processes in a body of water. 
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Physical processes include settling of sediments and substances, as well as phenomena such as aeration. 
Biochemical processes transform one substance into another through adsorption, decay, or biological 
processes. Sources add or subtract mass from the computational cell. Additional pollutant mass can come 
due to sewer wasteloads and subtraction is due to intakes. 

The three change processes described above (TR, P and S) are governed by the extended transport 
equations (Somlyódy & van Straten, 1986), which is written in three dimensions as: 


C əf, a) 9 à(, ac) a a(. CY 9 
a aE) a Oo [o a o S (P. | ag ©) 


+ rc(C, param) (8.3) 


where C(c,, C>,..., c,) is the n dimensional mass concentration vector for n constituents of analysed water 
quality; u(u,, u, u.) is the velocity vector; D(Dx, Dy, Dz) is a vector containing diffusion coefficients for 
the x, y and z, directions respectively; and rc(C, param) is the n dimensional vector of rates of change of 
constituents due to P and S processes as a function of concentrations of the n constituents and a number of 
parameters, param that depend on the model choice. Model calibration is usually performed to determine 
the values of param. 

In case sediment is a component of the water quality model a mass conservation equation for the 
sediment is added to the model along with corresponding interface terms sediment-water and water-air. 

Equation (8.3) is a system of PDEs that are solved numerically (together with the sediment equation in 
case sediment is considered). Another approach to solve the formed system of PDEs is to use a number of 
m interconnected conceptual elements (so-called tanks) that lead to a system of n x m ODEs. 

Models that are defined in steady state are the ones for which there are no changes in the concentrations 
with time, that is, the term 0C/dt in equation (8.3) equals zero. 

Temperature affects water quality processes, hence it is important to incorporate temperature and its 
modelling when it varies over the simulation period, or when the inflow to the water body has a heat that 
is different from the one of the water body that is modeled. Temperature models are based on heat balance 
of the modeled water body, by taking into account solar and atmospheric radiation and evaporation. Such 
models are not easy to set-up, because they require a lot of measured data that is not always available. 


8.3 RIVER WATER QUALITY MODELS 


River waters are usually receiving pollutants from several sources such as nutrients carried by sediments 
discharged in a river by runoff; sewers and storm sewer overflows; and so on. Mathematical models are 
useful in describing the state of water quality in a river system, as well as to predict the change in state, in 
case initial and/or boundary conditions are changed. Examples of changes are: changes in the type (point 
or non-point) of loads; morphological changes of the river bottom, due to construction and/or operation of 
control structures (weirs, dams, etc.); or changes in hydrograph inputs (upstream and lateral) due to effect 
of climate change. The degree of complexity in describing river water quality state varies from model to 
model. Initially, river water quality models determined dissolved oxygen concentrations (DO) and, later 
were extended to predict plankton-nutrient dynamics in river eutrophication problems (Rauch et al., 1998). 
Other constituents that may be incorporated in water quality models are macrophytes (Barendregt & Bio, 
2003; Zalewski, 2004); phosphorus; oil spills and pesticides. A comprehensive overview of existing river 
water quality models is given by Cox (2003). 
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Oxygen concentration of river waters is the most important factors affecting aquatic life. Hence water 
quality standards are related to the concentration of dissolved oxygen (DO) and parameters affecting 
it. Discharging pollutants in receiving river waters has as effect a reduction in the oxygen content. A 
related element is the Biochemical Oxygen Demand (BOD), which shows the oxygen consumption due to 
degradation of organic material. BOD5 is the amount of BOD consumed over 5 days. 

Oxygen consumption is described in form of a first order ordinary differential equation: 


dCgop 


dt e -(K, +K,s) i Crop t K(C,, m C,) = D, (8.4) 


where CBOD is the concentration of BOD (M/L3); K1 the decay constant (sec); KIS the sedimentation/ 
adsorption constant (sec); Cos the oxygen saturation constant (M/L?); C, the oxygen concentration in the 
river (M/L?); K the re-aeration kinetic constant (sec) and D, the rate of oxygen removal caused by benthic 
processes. Consequently oxygen concentration in a river is determined using the equation: 


dC 
His = —K, - Cgop (8.5) 
An important phenomena in modelling water quality in rivers is re-aeration that occurs at the air-water 
interface: 
IR 
mi x K(C,, = C,) (8.6) 


Removal of organic material due to sedimentation process (Dobbins, 1964) is represented in equation 
(8.4) by the sedimentation/adsorption constant. The actual amount of BOD sedimented/adsorbed at the 
bottom of a river was described by Harremoes in 1982, which is an important component of river water 
quality analysis when storm sewer overflows are involved because of their high content in suspended 
organic content ready to be settled. 

Oxygen concentrations are important for the aquatic life, however other factors are also important, 
such as is the case of ammonia in non-ionized form (NH?), which is toxic to fish. The pH value in a 
river determines the dissociation of ammonia into ionized and non-ionized, the higher the PH the higher 
fraction of ammonia is non-ionized. At pH values above 9 almost all ammonia is in a non-ionized form. 
Values of pH should not be too high, neither too low, because high pH shows increasing acidification, 
while low pH is toxic to the aquatic life, because alluminium’s solubility increases. 


8.4 LAKES WATER QUALITY MODELLING 


The above presented water quality modelling concepts are applicable as well to lakes, coastal or ocean 
waters. This section presents just few aspects that are specific to water quality modelling in lakes and 
reservoirs. Usually lakes and reservoirs are studied as water bodies with the same physical properties, 
despite their minor differences. Reservoirs are a result of a built structure (most often a dam), and are used 
primarily to address problems of water shortages, excess floods, or energy production (ILEC, 2005). 

If a river flows into a lake there are significant changes in the state of its water quality because its velocity 
reduces. All materials (pollutants and sediments) are carried at higher velocities in rivers and they enter 
to slow moving waters. In lakes and or reservoirs there are more opportunities for algae (phytoplankton) 
growth and development of eutrophication. 
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Water quality in lakes and reservoirs is assessed based on several characteristics of the water body. 
These were defined by Louks and van Beek (1983) as: 


* clarity or transparency of the water, that is, greater water clarity indicates better water quality; 
* concentration of nutrients, that is, lower concentrations show better water quality; 

* quantity of algae, that is, lower levels of algae indicate better water quality; 

* oxygen concentration, that is, the higher the concentration of oxygen the better for fish; 

* concentration of dissolved minerals, that is, low values give better water quality; 

* pH value (a neutral pH of 7 is advised). 


Supply of water to a lake is important for determining its water quality. There are two main sources of 
water supply for lakes and reservoirs; precipitation and river inflow. If precipitation is the major source 
of water, the lake is acidic and low in nutrients. If river inflows are the major source for water supply, lake 
water quality is more variable, and has high nutrients level. For the later case water quality depends on the 
volume of inflow and runoff from the surrounding slopes. Human activities around a lake have important 
effect on its water quality. 

Lakes and reservoirs have three fundamental characteristics that make them quite different in physical 
behaviour as compared with rivers. These characteristics include: mixing, long retention time, and response 
dynamics. 

Reservoirs or lakes receive inputs from diverse sources; precipitation, river inflow, heat and 
wind-induced energies that cause waves, thermal energies, contaminants, nutrients and organic substances; 
which are mixed within a lake so that both the resources and problems are dissipated throughout its 
volume. 

Lakes and reservoirs are generally slow to respond to changes due to their long retention times. They 
store floods, pollutants, and heat without experiencing immediate changes. The implications of the 
long retention times is that lakes are relatively stable; allow for suspended materials to settle to the 
bottom, acting as high efficiency sinks for many materials; and allow for complex, and sometimes unique 
ecosystems to evolve. The drawback of long retention time is that, once a lake or reservoir is degraded, it 
takes a long time for it to recover or be restored. 

The response dynamics of lakes to a nutrient concentration is different from one of a river. Lakes and 
reservoirs respond to changes in a highly non-linear way, as depicted in Figure 8.2. Degradation in response 
to increased nutrient concentrations (from A to B), may not be visible until nutrient concentrations are 
high and the water body abruptly switches its status. In Figure 8.2 the plankton concentrations are high, 
and the lake goes from status B to C. The lake cannot simply go back from C to B. There are likely to be 
irreversible changes to the ecosystem, so recovery follows a path from C down to A through D. 


» 


Plankton concentration 


Nutrient concentration 


Figure 8.2 Response of a lake to eutrophication process. 


Downloaded from https://iwaponline.com/ebooks/book-pdf/650788/wio9781780400457.pdf 
bv IWA Publishina. publications@iwap.co.uk 


154 Computational Hydraulics 


Water quality of a lake is affected by the extent to which the water mixes. The depth, size and shape 
of a lake are the most important factors influencing mixing, along with other factors such as climate, 
lakeshore topography and inflow from rivers. Density of water has highest value at 4?C. Variations 
in density caused by difference in temperatures prevent warm and cold water from mixing. Uniform 
water density facilitates complete mix of lake waters, recharging the bottom water with oxygen and 
bringing nutrients to the surface. Wind and waves circulates warm water only 6—8 m deep. If a lake 
is shallow the water may stay completely mixed for the entire summer season. During summer season 
deep lakes experience stratification into three zones: epilimnion (warm surface layer), thermocline or 
metalimnion (transition zone between warm and cold water), and hypolimnion (cold bottom water). 
Stratification traps nutrients from bottom sediments in the hypolimnion. In the fall season, the lake 
surface cools and eventually water temperature becomes uniform from top to bottom, hence mixing 
appears (fall overturn). During fall season often algae bloom appears to the surface of a lake because 
of the mixing. 


8.5 EXAMPLES OF LAKE HYDRODYNAMICS AND WATER 
QUALITY MODELS 


8.5.1 Sontea-Fortuna wetland system 


The Danube Delta is the largest wetland of the Danube river system (see Figure 8.3). It consists of 
a network of river channels, lakes, and wetlands that help to improve water quality of the receiving 
waters of the Danube, the Black Sea. Danube Delta forms a buffer zone that filters pollutants from the 
Danube river, while helping to reduce flood peaks by storing water. A management plan of the Danube 
Delta is in place in order to ensure its protection and to establish the preservation of an ecological 
equilibrium. Wetlands of the Danube Delta are part of this management plan, however their functioning 
is difficult to be quantified because of the dynamic nature of the system (Popescu et al., 2014). Status of 
functioning of a wetland in the Danube Delta is evaluated based on a clear defined number of indicators 
(EU WFD, 2000). These indicators are used to assess the condition of the ecosystem, by monitoring 
trends in time and space. Sontea-Fortuna wetland is one of the most important wetlands of the Danube 
Delta system because of its location at the entry point in the delta (see Figure 8.4). Terrestrial and 
aquatic habitat of this wetland depends on the water level in the Sontea-Fortuna area. The major 
driving force for changes in the system is the hydrological regime of the canal system of the wetland, 
including periodic floods. 

The Sontea-Fortuna wetland is located on Romanian part of Danube Delta and covers a surface of 
24,636 ha between Chilia branch to the West and North; Sireasa and Stipoc canals to the North; Catavaia 
canal and Old Danube to East; and Tulcea and Sulina branches to the South (see Figure 8.3). The area 
elevation of the wetland varies from as low as —6 m m.b.s.l to +4.5 m m.b.s.l. The lowest elevation in the 
area occurs along Sontea canal as the central axis of the Sontea-Fortuna system. 

Results of a 2D structured grid hydrodynamic model are shown here. The model is used is to determine 
the behaviour of the system during high and low flow conditions, in order to predict the hydrological 
regime in the Sontea-Fortuna wetland. The aim of the study is to determine the area exposed to dry 
conditions, especially during low flows. Three indicators shows the status of the Sontea Fortuna wetland; the 
inundation pattern, which is in a good status if the water level in the wetland ranges between 1.5 m m.b.s.1 
and 3.5 m m.b.s.1; flooding duration, that should not exceed 30 days and circulation of water in Fortuna 
lake, which should not be in a stagnant state. 
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Figure 8.4 Water levels in Sontea Fortuna wetland for dry hydrological conditions. 
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Model is setup for three hydrological conditions; wet, normal and dry year, which correspond to high, 
medium, and low water level respectively, based on measured data kindly provided through the EU FP7 
research project EnviroGRIDS (Popescu et al., 2014). The Sontea—Fortuna area is modelled using Delft3D, 
which is a tool developed by Deltares (Deltares, 2010) to simulate hydraulic phenomena in river, estuarine 
and coastal areas. The software simulates variations in time and space (2D or 3D) of hydrodynamics, 
morphology water quality and sediment transport phenomena. The Sontea—Fortuna grid is a combination 
of grid size cells that varies from 90 m x 60 m to 160 m x 90 m. Total number of computational cells is 
37039. 

Figure 8.4 shows water levels in the area, as results of the model, for dry hydrological conditions. 
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Figure 8.5 Water levels at upstream boundary of Sontea-Fortuna wetland for dry hydrological conditions. 


Results help in identifying dry areas, where habitat is at risk and water quality conditions deteriorate, for 
too low water levels. Treshold for water levels are set by policies of the Romanian Ministry of Environment. 

Figures 8.5, 8.6 and 8.7 shows results of water levels predictions in three points of the Sontea-Fortuna 
wetland for the three modelling cases; wet condition (October-December, 2010), dry condition (June— 
August, 2010) and one adaptation measure for the case of dry condition. The three points are located in the 
upstream, middle and downstream of the wetland, respectively. 

For the three considered cases it can be seen that water level can be maintained between and above the 
levels of 1.5 m (m.b.s.1) to 3.5 m (m.b.s.I), and flooding does not exceed 30 days. 

Based on the outcomes of the model adaptation measures can be proposed for improving the ecological 
status of the wetland. 


8.5.2 Lake Taihu water quality 


Lake Taihu is a shallow, freshwater wind-driven lake in the Yangtze Delta, located at the border between 
the provinces of Jiangsu and Zhejiang in Eastern China (see Figure 8.8). Eutrophication is always a serious 
problem in Taihu Lake. Water pollution associated with population increase, massive industrial growth 
with highly polluting heavy industry has taken a serious toll throughout China (Boqiang, 2007). 
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Figure 8.6 Water levels in the middle of Sontea-Fortuna wetland for dry hydrological conditions. 
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Figure 8.7 Water levels at downstream boundary of Sontea-Fortuna wetland for dry hydrological conditions. 
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Figure 8.8 Lake Taihu location. 
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Many studies have been carried out regarding the water quality condition in the Taihu Basin. Zhu and 
Geng (2004) show that according with Chinese surface water quality standards there are 5.99 billion n? 
of water (4296 in the total runoff) that are severely polluted and unusable for any user in Lake Taihu. 
Wang et al. (2011) estimated the pollutant fluxes of rivers in Zhejiang Province as part of the Taihu Lake 
Basin and indicated that the increasing levels of inflow of pollutants in the lake is the main reason for the 
deterioration of water quality in the lake. When comparing the mean values to those obtained between 1991 
and 2001, the water quality parameters such as TN and TP have had higher values in the recent five years. 
Mao et al. (2008) published a three-dimensional eutrophication model of Lake Taihu, which fully couples 
the biological processes and hydrodynamics, and takes into account the effects of sediment release and the 
external loads from the tributaries. The main outcome of such a model is that two of the Lake Taihu's bays; 
Zhushan Bay and Xu Bay; are susceptible to algal blooms with high chlorophyll-a concentration. 
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Figure 8.9 Spatial distribution of Nitrogen and Phosphorus in Taihu lake, during year 2008. 
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Here are presented water quality model results using a two-dimensional hydrodynamic model coupled 
with a two-dimensional water quality model. Model is set-up to investigate possible effects of wind on 
the eutrophication of the lake. A Delft3D model is built to simulate the hydrodynamic state of the lake. 
Hydrological, topographical and water quality data to set-up an example model are collected from published 
papers of Yue et al. (2013); Hu et al. (1998); Zhu and Geng (2004) and Yiping and Zhougho (2011). Three 
different sets of wind conditions (no wind, constant wind and measured wind) are considered as wind 
effects on the water level and velocity, in wet and dry seasons. Based on the calibrated hydrodynamic 
model a two dimensional water quality model is built and calibrated based on estimated data inputs for 
nutrient loads. A one-year hydrodynamic simulation is performed, for the year 2008, in order to show 
possibilities for determining the state of water quality in Lake Taihu. The estimated nutrient loads included 
total nitrogen (TN) and total phosphorus (TP) concentrations. The estimates are based on the population 
of cities around Lake Taihu, and the wastewater treatment plant capacity of the provinces of Jiangsu and 
Zhejiang, where the main sources of pollutants come from. Two time periods (spring and autumn) of algae 
blooms development for the year 2008 are noticed to form in the lake. 

Figures 8.9 and 8.10 shows different modelling results for nitrogen, phosphorus and algae bloom for 
Taihu Lake, for the year 2008, for conditions of constant wind conditions. 
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Figure 8.10 Temporal distribution of Algae bloom in Taihu lake, during year 2008. 


The results of the study indicate that the hydrodynamic model with constant wind conditions (southeast, 
5 m/s) shows better results than the one where measured wind conditions are used, simply because the 
measured wind data is insufficient; and an increase in population in the area would lead to an increase in 
the maximum value of algae concentration from 35% to 76%. An improvement of technologies used for the 
treatment of the wastewater may bring the maximum value of algae concentration between 17% and 42%. 
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m E m T m 


Computational Hydraulics introduces the concept of modeling and 
the contribution of numerical methods and numerical analysis to 
modeling. It provides a concise and comprehensive description 
of the basic hydraulic principles, and the problems addressed 
by these principles in the aquatic environment. Flow equations, 
analytical and numerical solutions are included. 


The necessary steps for building and applying numerical methods 
in hydraulics comprise the core of the book and this is followed by 
two different example applications of computational hydraulics: 
river systems and water quality modelling of lakes and rivers. 


The theory and exercises included in the book promote learning 


of concepts within academic environments. 


Computational Hydraulics is intended for under-graduate and 
graduate students, researchers, members of governmental 
and non-governmental agencies and professionals involved in 


management of the water related problems. 
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